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ABSTRACT 

A new generation of ground-based interferometric detectors for gravitational waves is currently under 
construction or has entered the commissioning phase (LIGO, VIRGO, GEO600, TAMA). The purpose 
of these detectors is to observe gravitational waves from astrophysical sources and help improve our 
understanding of the source origin and physical properties. In this paper we study the most promising 
candidate sources for these detectors: inspiraling double compact objects. We use population synthesis 
methods to calculate the properties and coalescence rates of compact object binaries: double neutron 
stars, black hole-neutron star systems and double black holes. We also examine the formation channels 
available to double compact object binaries. We explicitly account for the evolution of low-mass helium 
stars and investigate the possibility of common-envelope evolution involving helium stars as well as two 
evolved stars. As a result we identify a significant number of new formation channels for double neutron 
stars, in particular, leading to populations with very distinct properties. We discuss the theoretical and 
observational implications of such populations, but we also note the need for hydrodynamical calculations 
to settle the question of whether such common-envelope evolution is possible. We also present and discuss 
the physical properties of compact object binaries and identify a number of robust, qualitative features as 
well as their origin. Using the calculated coalescence rates we compare our results to earlier studies and 
derive expected detection rates for the Laser Interferometer Gravitational- wave Observatory (LIGO). 
We find that our most optimistic estimate for the first LIGO detectors reach a couple of events per year 
and our most pessimistic estimate for advanced LIGO detectors exceed ~ 10 events per year. 

Subject headings: binaries: close — gravitational waves — stars: evolution, formation, neutron, black 
hole 



1. INTRODUCTION 

Gravitational waves are a natural consequence of Ein- 
stein's theory of general relativity (Einstein 1916, 1918). 
Indirect evidence for their existence came first from ob- 
servations of the orbital decay of the Hulse- Taylor binary 
pulsar (Hulse & Taylor 1974, 1975a, 1975b; Taylor & Weis- 
berg 1982, 1989). Direct detection though and analysis 
of gravitational-wave sources are expected to provide a 
unique insight to one of the least understood of the fun- 
damental forces. They will also allow us to investigate the 
physical properties of astronomical objects that so far have 
been elusive because they do not emit any electromagnetic 
radiation (e.g., double black holes in isolation). 

A number of interferometers designed for gravitational- 
wave detection are currently in operation, being devel- 
oped, or planned. First results from the Japanese instru- 
ment TAMA300 have already been reported (Tagoshi et 
al. 2001). More advanced ground-based observatories are 
under development or in the commissioning phase, i.e., 
the US Laser Interferometer Gravitational-wave Observa- 
tory (LIGO; Abramovici et al. 1992), the British/German 
GEO-600 Observatory (Danzmann et al. 1995), the 
French/Italian VIRGO project (Caron et al. 1995), and 
the Australian ACIGA project (Sandeman 1998). Also 
future space missions are being planned, such as Laser In- 
terferometer Space Antenna (LISA; Bender et al. 2000), a 



joint project of ESA and NASA with the launch planned 
for 2009. 

Astrophysical sources of gravitational radiation relevant 
to ground-based interferometers include: inspiraling dou- 
ble compact objects, binary stars, rotating neutron stars, 
neutron star instabilities, supernovae, supermassive black 
holes and stochastic background (for a review see Thorne 
1987). Some of the most promising candidates are the in- 
spiral and coalescence of double compact objects (DCO), 
such as NS-NS, BH-NS, and BH-BH binaries. Success- 
ful detection of these sources at reasonable event rates 
depends not only on the instrument sensitivity and the 
strength of the gravitational- wave signals, but also on the 
coalescence rates and the physical properties of the sources 
out to the maximum distances of reach of a given instru- 
ment. 

Coalescence rates of double compact objects have been 
investigated by a number of different groups. For NS- 
NS binaries the rates can be calculated in two ways: (1) 
theoretically, based on the predictions of binary popula- 
tion synthesis calculations (e.g., Bethe & Brown 1998; 
Portegies Zwart & Yungelson 1998; Bloom, Sigurdsson 
& Pols 1999; Belczynski & Buhk 1999; Fryer, Woosley 
& Hartmann 1999) or (2) empirically, based on the ob- 
served sample of Galactic binary pulsars (e.g., Narayan, 
Piran, & Shcmi 1991; Phinney 1991; Curran & Lorimer 
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1995; Arzoumanian, Cordes & Wasserman 1998; Kalogera 
et al. 2001; Kim, Kalogera, & Lorimer 2001). At present 
both methods appear burdened with significant uncertain- 
ties, the theoretical approach due to the many poorly con- 
strained evolutionary model parameters, and the empirical 
estimates due to small-number sample of observed NS-NS 
systems (for a more detailed discussion see Kalogera et al. 
2001). Predictions for BH-NS and BH-BH binaries can be 
obtained only from theoretical calculations, since no such 
systems have yet been observed. Accurate predictions for 
BH binaries are even harder to obtain, as the evolution of 
high-mass stars still challenge our understanding. Several 
results of population synthesis calculations have been pre- 
sented in the literature, for BH-NS systems: Fryer et al. 
(1999); Belczynski, Buhk & Zbijewski (2000); Portegies 
Zwart & Yungelson (1998), and for BH-NS and BH-BH 
systems: Lipunov, Postnov & Prokhorov (1997); De Bon- 
der & Vanbeveren (1998); Nelemans, Yungelson & Porte- 
gies Zwart (2001b). 

Given the history of population synthesis calculations, 
one may wonder what is the purpose on another set of 
models of the formation of double compact objects. The 
primary goal of this present study is to investigate the 
relevant formation processes in the light of recent devel- 
opments in the evolution of compact object progenitors in 
binaries, and to clearly identify classes of formation chan- 
nels, their origin, and their relative contributions. We fo- 
cus not just on predicted rates, but at some level more 
importantly on the physical properties of formed binaries. 
Our parameter study is targeted to uncovering the sys- 
tematic uncertainties as well as identifying those (quali- 
tative) properties that appear to be robust. Investigation 
of these properties will help in characterizing the range 
of possible gravitational-wave signals and identification of 
their origin will allow reevaluating our results in the fu- 
ture as our understanding of model unknowns improves. 
As part of our analysis we also examine the implications 
of common-envelope phases involving helium stars and as 
a result we identify a number of new formation channels 
with important implications for the properties of double 
compact object^. Even though stellar structural charac- 
teristics may tentatively support such common-envelope 
phases late in the evolution of helium stars, we note that 
detailed hydrodynamical calculations will be needed to as- 
sess the viability of such phases involving helium stars or 
pairs of evolved stars as it has been assumed also in earlier 
studies (see Fryer et al. 1999; Nelemans et al. 2001a) 

Our paper is organized as follows: In § 2, we present 
a detailed description of our population synthesis code 
(named StarTrack). In §3 we: describe the extent of our 
parameter study (§3.1), discuss the statistical accuracy 
of our calculations (§3.2), present the formation paths of 
double compact objects (§ 3.3), and analyze in detail their 
properties (§ 3.4) and coalescence rates (§ 3.5) for a large 
number of different models. In § 4, we discuss the impli- 



cations of the most important of our results and conclude 
with prospects for gravitational-wave detection by ground- 
based interferometers. 

2. POPULATION SYNTHESIS MODEL 

2.1. Single- Star Evolution 

2.1.1. Overview 

In all the population synthesis models, we employ the 
analytic formulae derived by Hurley, Pols & Tout (2000, 
hereafter HPT) to model the evolution of single stars. 
These formulae represent fits to results from a large num- 
ber of stellar evolution models calculated by Pols et al. 
(1998), using the Eggleton stellar evolution code (Eggle- 
ton 1971,1972,1973; latest updates and input physics de- 
scribed in Han, Podsiadlowski & Eggleton 1994 and Pols 
et al. 1995). 

The HPT formulae^ allow us to calculate the evolu- 
tion of a star given its Zero Age Main Sequence (ZAMS) 
mass and its metallicity {Z). Stellar evolution is followed 
from ZAMS through a sequence of evolutionary phases 
depending on the initial (ZAMS) stellar mass: Main Se- 
quence (MS), Hertzsprung Gap (HG), Red Giant Branch 
(RG), Core Helium Burning (CHeB), Asymptotic Giant 
Branch (AGB), which is further divided into Early AGB 
(EAGB) and Thermally Pulsing AGB (TPAGB), and for 
stars stripped off their hydrogen-rich layers: Helium Main 
Sequence (HMS), Hehum Giant Branch (HGB). We end 
the evolutionary calculations at the formation of a stellar 
remnant: a white dwarf (WD), a neutron star (NS) or a 
black hole (BH). During the evolution we can track some 
of the basic stellar parameters: radius, luminosity, stellar 
mass, and core mass. 

These single star models include the effects of mass loss 
due to stellar winds as described in HPT. Mass loss rates 
are adopted from the literature for different evolutionary 
phases. For H-rich stars on MS (Nieuwenhuijzen & de 
Jager (1990); using Z dependence of Kudritzki et al. 1989), 
for RGB (Kudritzki & Reimers 1978) and AGB (Vassil- 
iadis & Wood 1993) stars, and for Luminous Blue Vari- 
ables (HPT). For He-rich stars W-R, mass loss is included 
using rates derived by Hamann, Koesterke & Wessolowski 
(1995) and modified by HPT. Given the importance of sin- 
gle star winds in the formation of compact objects (e.g.. 
Brown, Weingartner, & Wijers 1996; Ergma & van den 
Heuvel 1998; Fryer & Kalogera 2001), we examine this el- 
ement of single star evolution and calculate synthesis mod- 
els with varying strengths of stellar winds. 

We have introduced two modifications to the HPT for- 
mulae concerning the treatment of (i) Helium-star evolu- 
tion, and (ii) final remnant masses. These are described 
in detail in the next two subsections. 

2.1.2. Helium Star Evolution 



^After the submission of this paper to ApJ, an early paper by Tutukov & Yungelson (1993) came to our attention. In this paper the authors 
allow for mass transfer from helium stars without a merger, although not through CE evolution. There is qualitative agreement between their 
and our results, however we point out that there are many differences in the binary evolution modeling (e.g., Tutukov & Yungelson do not 
account for any kicks to NS or BH not for hypercritical accretion, and many more). We also found out that CE evolution for low-mass helium 
stars has been discussed qualitatively in the past, in the context of the formation of NS-NS systems and X-ray binaries (van den Heuvel 1992; 
Taam 1996). 

^We implemented the formulae into a C code, compared the results of our implementation with the original HPT Fortran subroutines for 
ZAMS masses 0.5— 100 Mq and for metallicities: Z = 0.0001 — 0.03, and found perfect agreement. 
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After core helium exhaustion low-mass helium stars ex- 
pand significantly and develop a "giant-like" structure 
with a clearly defined core and a convective envelope (Del- 
gado, & Thomas 1981; Habets 1987; Avila-Reese 1993; 
Woosley, Langer, & Weaver 1995). Radial expansion is 
very important for stars in binaries as it may lead to mass 
transfer episodes. Low-mass evolved helium stars with 
convective envelopes possibly transfer mass on a dynamical 
time scale, and as a consequence a common envelope (CE) 
phase may ensue^. On the other hand, relatively massive 
helium stars do not develop deep convective envelopes and 
do not experience any significant radial expansion. Even if 
they happen to initiate mass transfer (MT), they respond 
differently to mass loss and a CE phase is not likely. 

Based on the HPT formulae, the implied upper mass 
limit for helium stars to develop deep convective envelopes 
is Mconv ~ 2.2 Mq), which we consider to be rather low. 
We have examined in detail models of evolved helium stars 
(Woosley 1997, private communication) and found that he- 
lium stars below 4.0 Mq have deep convective envelopes, 
whereas slightly more massive helium stars (~ 4-4.5 Mq) 
still form convective envelopes, although shallower. As an 
example, an evolved helium star of M = 2.5 Mq devel- 
ops a deep convective envelope reaching down almost to 
the core through 90% of the stellar radius. Based on this, 
for our standard model, we adopt Mconv = 4.5 Mq for an 
evolved helium star to develop a deep convective envelope. 
However, in face of the uncertainties and the importance 
of helium-star evolution, we treat this value as a model pa- 
rameter and include it in our parameter study. Through- 
out this paper we refer to the evolved helium stars with 
M < Mconv as low-mass helium stars/giants (or low-mass 
HGB stars) and to the stars with M > Mconv as massive 
helium stars/giants (or massive HGB stars). We note that, 
apart from the issue of convective envelopes, a compari- 
son in terms of radial evolution alone between the HPT 
formulae and other published models of helium stars (e.g. 
Habets 1987) shows excellent agreement, so we adopt these 
formulae in calculating stellar radii. Detailed discussion of 
helium star evolution in context of double compact objects 
formation is given in Belczynski & Kalogera (2001). 

2.1.3. Remnant Masses 

The masses of NS and BH calculated by HPT are very 
small. For progenitors of M^ams = 40 — 100 Mq, they ob- 
tain final remnant masses in the range ~ 1.8 — 2.0 Mq, 
(see their Figure 20). On the other hand, measurements 
of BH masses in binaries, although still highly uncertain, 
yield much higher values: ^ 3-20 Mq (Orosz et al. 2001; 
McClintock et al. 2001; Froning & Robinson 2001; Wag- 
ner et al. 2001; Table 1 of Fryer & Kalogera 2001 and 
references therein). The HPT analytic formulae are based 
on the stellar models evolved at least to the formation of 
a carbon-oxygen (CO) core, and these core masses are in 
good agreement with earlier stellar models (e.g., Schaller 
et al. 1992). Therefore, instead of using their somewhat 
arbitrary prescription (see their equation [75]) for NS and 

^We note that whether a common-envelope phase develops as well as the outcome depends not only on the presence of convective envelopes 
but also on the detailed stellar structure (density profile, masses of the envelope relative to the core) as well as the mass ratio in the binary 
and the spin evolution of the envelope during the mass-transfer episode (Fryer 2001; Rasio 2001; Taam 2001, private communications). The 
hydrodynamics of mass transfer from helium stars with appropriate masses and companions have never been studied so far, but such studies 
will be necessary to determine the outcome of these evolutionary phases. In the present study we account for such a possibility and examine 
its implication for the DCO populations. 



BH masses, we adopt a different formula derived on the 
basis of core-collapse hydrodynamical calculations (Fryer 
1999, Fryer & Kalogera 2001). We note that there are stiU 
many uncertainties associated with the physics of the su- 
pernova mechanism, and there is still disagreement among 
research groups about the cause of the explosion (e.g., 
Liebendoerfer et al. 2001 and references therein). For our 
study we choose the specific core-collapse models, primar- 
ily because they are the only one available in the literature 
that allow us to calculate the remnant mass taking into ac- 
count fallback, in a self-consistent way. 

We adopt the HPT formulae to calculate the final CO 
core mass of a star, and we use stellar models of Woosley 
(1986) to obtain a final FeNi core mass as a function of 
the CO core mass. At the time of core collapse, we de- 
scribe the star with its current mass M, the final CO core 
mass Mco and the final FeNi core mass MpoNi- To esti- 
mate the mass of the remnant formed, we follow the re- 
sults of hydrodynamical calculations of core collapse (Fryer 
1999). These calculations have shown that (i) progenitors 
with Mzams < 20 Mq do not experience any significant fall 
back and apart from the collapsing core, the stellar outer 
layers are expelled as SN ejecta, (ii) for progenitors with 
Mzams > 42 Mq, the whole pre-coUapse star directly im- 
plodes to form a BH and no SN event takes place, (iii) 
for intermediate progenitor masses remnants are formed 
through an initial core collapse and subsequent partial fall 
back. The outcome of a given core-collapse event is found 
to depend primarily on the core mass of the collapsing star 
and not on the initial (ZAMS) stellar mass. Taking into 
account that the stellar population models used by Fryer 
(1999) were of constant mass (without wind or mass trans- 
fer effects; see Fryer & Kalogera 2001), we have developed 
a prescription based on CO core masses. Using the HPT 
formulae we obtain Mco = 5Mq for M^ams = 20 Mq, and 
Mco = 7.6 Mq for M^ams < 42 Mq, for Z = 0.02 using the 



standard wind mass loss prescription, 
the mass of the remnant as follows: 



Mr, 



MpoNi 

MpoNi + /fb(M - 

M 



MpcNi) 



We then calculate 



Mco < 5 Mq 
5 < Mco < 7.6 
Mco > 7.6Mq 



(1) 

where /fb is the fall-back factor, i.e. the fraction (from to 
1) of the stellar envelope that falls back. In Figure |l| we 
present the initial-remnant mass relation, for single Pop- 
ulation I stars {Z = 0.02), and for three different wind 
mass-loss rates. 

The top panel of Figure |l| shows the initial-remnant 
mass relation for our standard prescription of wind mass- 
loss rates. 

Stars of Afzams — 0.5-0.8 Mq form hehum white dwarfs 
(He WD), of Mzams = 0.8-6.3 Mq form carbon-oxygen 
white dwarfs (CO WD), and of M^ams = 6.3-8.0 Mq form 
oxygen-neon white dwarfs (ONe WD). 

For stars with initial mass M^ams in excess of 8 Mq, 
the stellar core mass exceeds the Chandrasekhar mass and 
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once nuclear reactions stop the core collapses into a NS. 
For Mzams = 8 — 20.0 Mq, the outer stellar layers (on av- 
erage ~ 7Mq) are lost in SN explosion, and no fall back 
takes place. This part of the initial-remnant mass relation 
is characterized by a slow increase of the NS mass with 
increasing Mzams, as more massive progenitors form only 
slightly more massive cores, which eventually collapse to 
NS. 

Partial fall back is initiated above Afzams — 20.0 
and the amount of material accreted back onto the collaps- 
ing core increases linearly with increasing progenitor mass. 
Due to the fall back the initial-remnant mass relation sig- 
nificantly steepens over Mzams = 20.0 Mq. Note that the 
remnant mass reaches maximum NS mass M^fx — 3.0 Mq 
at Mzams = 20.7 Mq, SO in this model the most massive 
NS are formed through fall back. 

For Mzams ~ 25 Mq there is a small drop and a subse- 
quent flattening of the initial-remnant mass relation. At 
this and higher progenitor masses stellar winds strip star 
completely of its hydrogen layers. Once a naked helium 
star is formed the wind mass loss rate increases signifi- 
cantly. The effect of wind mass loss enhancement with in- 
creasing Mzams diminishes both the final core mass and the 
pre-coUapse envelope mass (the fall back mass reservoir) 
and therefore decreases the remnant mass. With increas- 
ing Mzams the fall back factor /fb goes up, but because of 
the helium star mass loss rate effect on the pre-collapse 
star, the remnant mass does not increase as fast as before. 

Stars more massive than M^ams = 42 Mq collapse di- 
rectly into a BH at the end of their evolution. The initial- 
remnant mass relation flattens out even more. This is a 
consequence of a "fall-back saturation" . The fall back is 
now complete (/fb = 1) and cannot increase any more with 
progenitor mass, as it was the case for less massive stars 
in the range of partial fall back. Now the remnant mass 
depends only on the pre-collapse mass, which increases 
rather slowly with progenitor mass as the wind mass loss 
rate increases also with A/^ams- 

Stars over M^ams ^ 48 Mq reach high luminosi- 
ties and radii exceeding the Humphreys-Davidson limit 
(Humphreys & Davidson 1994) and enter the Luminous 
Blue Variable (LBV) phase. During the LBV phase stars 
lose mass at a very high rate, which result in a sudden drop 
in the initial-remnant mass relation. Due to this signifi- 
cant mass loss, stars over M^ams = 49 Mq form smaller 
cores and experience partial fall back at core collapse. 
The effect of LBV-like mass loss takes over the effect of 
an overall increase of CO core mass (and thus remnant 
mass) with Mzams, and it influences the remnant mass in 
two ways. First, LBV-like mass-loss rates increase with 
increasing stellar mass leading to a decrease of remnant 
masses. Second, for higher stellar ZAMS masses, the LBV 
phase begins earlier in the life of a star. Hence, the stars 
are stripped of their hydrogen layers earlier, and less mas- 
sive helium stars are formed, leading to a decrease in rem- 
nant masses. In the standard wind model, remnant masses 
reach their local minimum at M^ams = 52.5 Mq, when 
the LBV phase sets in the earliest possible, i.e., at the 
beginning of Hertzsprung gap. For more massive stars, 
the overall increase of remnant mass with progenitor mass 
dominates once again, as CO core masses steadily grow 
and the contribution of fall back steadily increases. 



For progenitors with Mzams > 72.3 Mq remnants are 
once again formed through complete fall back and the 
initial- remnant mass relation flattens out, following the 
slow increase of pre-collapse mass with Mzams- 

In the middle panel of Figure |], results are shown for 
Wolf-Rayet mass-loss rates reduced by a factor of two. 
The initial-remnant mass relation changes only slightly 
from that of the top panel as partial fall back at high 
ZAMS masses is eliminated. As expected, compact ob- 
jects are slightly more massive and the maximum BH mass 
increases from ~ 11 Mq (top panel) to ^ 15 Mq (see also 
Wellstein & Langer 1999). 

In the bottom panel of Figure ^ results are shown for 
the case of all wind mass-loss rates reduced by a factor of 
two, and we find some significant changes in the initial- 
remnant mass relation. First, as expected, the compact 
object masses are higher than for both models described 
above; the maximum mass of the compact object in this 
model increases to ~ 19 Mq. Second, the qualitative shape 
of the relation changes. As in the standard wind model, 
we find two sudden drops in the relation, with the first one 
being more prominent than before. In the standard model 
the first drop corresponds to the stripping star of its hydro- 
gen layers and the formation of a helium star, whereas the 
second drop is related to the onset of LBV phase for very 
luminous and extended stars. These two effects lead to sig- 
nificantly increased mass loss rates, and thus decrease the 
remnant masses. Detailed examination of this model re- 
veals that in the bottom panel, the first drop corresponds 
to the onset of the LBV phase and the second to the he- 
lium star formation. Since all wind rates are reduced and 
wind mass loss rates depend strongly on the stellar mass, 
helium stars are exposed only for stars of higher initial 
mass. Consequently, stars with much lower initial masses 
(below Mzams ~ 40 Mq) retain their hydrogen-rich layers 
and are able to evolve to higher luminosities and larger 
radii, which shifts the onset of the LBV phase to lower 
initial masses. Another difference is that, the immediate 
progenitors of compact objects tend to be more massive, 
and therefore direct BH formation spans much wider range 
of initial masses than for the standard wind model. These 
more massive progenitors form more massive CO cores, 
and none of the factors that decrease pre-collapse masses 
are strong enough to decrease the CO core masses below 
7.6 Mq (see eq. |l|). Consequently, BH formation through 
partial fall back does not occur for two separate Mzams 
ranges (as is the case for Mzams= 21-42 and 48-72 Mq in 
the standard wind model). 

2.2. Binary Evolution 

2.2.1. Overview 

We use Monte Carlo techniques to model the evolution- 
ary history of double compact objects and study their 
physical properties. We generate a large number {N > 
10^) of primordial binaries, drawing their initial phys- 
ical parameters (component masses, orbital separations 
and eccentricities) from assumed distributions described in 
§ 2.2.2. Besides that, Monte Carlo events are used to i) set 
the time of birth of each binary in the disk of our Galaxy 
(for constant star formation rate the probability distribu- 
tion of the birth time is flat in the range — lOGyr), ii) 
set the kick velocity direction and magnitude each NS or 
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BH receives when formed in an asymmetric SN explosion 
and in) set the position along the binary orbit in which SN 
take place. In all other instances, the evolution is modeled 
based on a set of formulae and prescriptions (described in 
what follows) that depend on the binary properties. 

The evolution of each system is followed in a number of 
time steps. The time step is chosen to be a small fraction 
(< 0.01) of the evolutionary phase lifetime of the more 
rapidly evolving stellar component, making sure that the 
relative change in radius is always less than 10%. We start 
the evolution of a given binary with two components at 
ZAMS, and we stop the calculation when the two compo- 
nents have formed stellar remnants or when a merger oc- 
curs or when the system has reached the present time, i.e., 
the sum of its age and its birth time is equal to lOGyr). 
The evolution of each binary component in isolation is 
calculated as described in § 2.1, and a number of evolution 
effects on the binary orbit (e.g., mass and angular momen- 
tum losses due to stellar winds) are taken into account. At 
every evolutionary time step, we check for possible binary 
interactions. If any of the components fills its Roche lobe, 
we calculate the effects of mass transfer and possible mass 
and angular- momentum loss on the orbit. Depending on 
the masses of the components and their evolutionary stage, 
we apply appropriate prescriptions for the modeling of dif- 
ferent mass transfer events (as described in § 2.2.4). If the 
binary survives the mass-transfer event (i.e., both stellar 
components fit within their Roche lobes), we continue to 
model its evolution. At every time step, we also check 
the single-star models (see §2.1) to determine if any of 
the components has finished its nuclear evolution, and has 
turned into a compact remnant. If a remnant is born in a 
SN explosion, we calculate the effects of SN kicks and mass 
loss on the binary orbit (both circular and eccentric). If a 
binary is disrupted in SN explosion, we follow the single- 
star evolution of components, otherwise we continue with 
the evolution of the whole binary. Once a binary con- 
sists of two remnants, wc; calculate its merger lifetime, i.e., 
the time until the components merge due to gravitational 
radiation and associated orbital decay, and we study the 
properties and formation rates of different classes of bi- 
naries containing compact objects. In what follows we 
describe in detail prescriptions for binary evolution and 
all model assumptions as chosen for our standard model. 
The assumptions for models in our parameter study are 
given in § 3.1. 

2.2.2. Distributions of Initial Parameters 

A binary system is described initially by four parame- 
ters: the mass of the primary Mi (the initially more mas- 
sive component), the mass ratio 9 = where M2 is the 
mass of the secondary (initially less massive component), 
the semi-major axis A of the orbit, and the orbital eccen- 
tricity e. We assume that the initial distributions of these 
parameters are independent. 

For both, single stars and binary system primaries, we 
adopt the initial mass function derived by Scalo (1986), 

*(M) oc M-^-\ (2) 

in the mass range Mi = 5 — 100 Mq relevant to compact 
object formation. Although, for single stars, the mini- 
mum initial mass for NS formation is ~ 8 Mq (or maybe 



higher), a lower limit of 5 Mq ensures that we do not miss 
any NS progenitors due to binary evolution. Mass trans- 
fer in a binary can increase any component's mass, and 
effectively decrease the minimum initial mass of NS pro- 
genitors. Let us consider a border-line case of the lowest- 
mass primary (SM©) and the corresponding highest-mass 
secondary (SM©). The primary evolves off the MS with a 
helium core of ~ 1 M© and ~ 4 Mq in its envelope. For our 
standard model of binary evolution, we assume that only 
half of the donor envelope can be accreted by its compan- 
ion (see §2.2.4). So for this border-line case, the initial 
secondary mass can be increased to < SM©, or for all 
other cases, to masses smaller than the mass limit for NS 
formation. 

Following Kuiper (1935), we assume a flat mass ratio 
distribution, 

$(g) = 1 (3) 

in the range g = — 1. Given value of the primary mass 
and the mass ratio, we obtain the mass of the secondary 
M2 = qMi. We evolve only systems with secondaries 
of Mzams > 0.5 Mq, as less massive stars will not evolve 
within a Hubble time. 

The distribution of the initial binary separations is as- 
sumed to be flat in the logarithm (Abt 1983), 

T{A) cx i, (4) 

where A ranges from a minimum value, such that the pri- 
mary fills its Roche lobe, up to 10^ R©. 

We adopt the thermal-equilibrium eccentricity distribu- 
tion for initial binaries, 

S(e) = 2e, (5) 

in the range e = — 1 (e.g., Heggie 1975; Duquennoy & 
Mayor 1991). 



2.2.3. Evolution of Binary Orbit 

Tidal circularization of binary orbits takes place in sys- 
tems in where the size of any component is comparable to 
the binary separation (Zahn 1978). This condition is sat- 
isfied when at least one of the binary components evolves 
beyond the MS or the initial binary separation is small. 
To calculate the effects of tidal circularization on a bi- 
nary orbit we follow the prescription developed by Porte- 
gies Zwart & Verbunt (1996): circularization occurs, if the 
stellar radius of one component is larger than 0.2 of the 
periastron binary separation. The orbital elements {A, e) 
change under conservation of angular momentum until the 
new periastron distance is equal to 5 stellar radii of the 
component driving the circularization, or until the orbit is 
circularized (e = 0). We also assume that tidal circulariza- 
tion takes place instantaneously (i.e. in one evolutionary 
time step). 

Stellar winds affect binary orbits through mass and 
angular- momentum losses. Assuming a spherically sym- 
metric stellar wind, which carries away the specific angu- 
lar momentum of the mass-losing star, and a circular orbit 
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(Jeans- mode mass loss), the change in the binary separa- 
tion is given by 



A{Mi + A/a) = const. 



(6) 



For eccentric orbits the change in A is rather similar (Van- 
beveren, Van Rensbergen & De Loore 1998). 

Orbits are also affected by mass transfer events and SN 
explosions. The treatment of these effects is described in 
the subsections that follow. 



2.2.4. Mass Transfer Events 

When one of the binary components fills its Roche lobe, 
mass is transferred to its companion through the inner 
Lagrangian point. The responses of both components and 
the orbit as well as the outcomes of MT phases depend on 
the masses and evolutionary stages of the two stars. In our 
calculations we distinguish between (i) dynamically stable, 
in general non-conservative (allows for mass and angular- 
momentum loss from the system) mass transfer, and (ii) 
dynamically unstable mass transfer that leads to common 
envelope evolution (for more details and reviews see, e.g., 
Ostriker 1975, Paczynski 1976, Iben & Livio 1993, Rasio 
& Livio 1996; Taam & Sandquist 2000). We are interested 
in the end products of MT events, and therefore, we as- 
sume that MT takes place instantaneously. Following any 
MT event, we update the binary component masses and 
evolutionary stages and calculate their radii ^. Using the 
post-MT radii and separation, we examine whether sys- 
tems survived the MT event and a component merger was 
avoided. 

Non- conservative Stable Mass Transfer. This phase is 
implemented if (i) the donor is a MS or HMS star trans- 
ferring mass to an accretor of any evolutionary stage, or 
(ii) the donor is a H-rich or He-rich giant and the accretor 
is not, and the following is true 



Mdon < c,M^ 



(7) 



where Mdon , -Mace ^rc the donor and accretor masses at 
the beginning of the MT episode. Based on earlier re- 
sults obtained by Hjellming & Webbink (1987), Kalogera 
& Webbink (1996) and Ritter (1999), Cr = 2.5 if the donor 
is in HG and the accretor is a MS star, and Cr = 1 in all 
other cases. 

During this type of MT episodes, part of the mass lost 
by donor (/a) is accreted onto the companion, and the rest 
is lost from the system with a specific angular momentum 
equal to (Podsiadlowski, Joss & Hsu 1992). The 

corresponding orbital separation A change can be calcu- 
lated for /a = from: 



A, 



2(J-1) 



Ml 



Ml 



M\ 



exp 



MLn) 



(8) 



and for /a > from: 

Af _Mi^ + Ml, 



M 



don 



acc 



Mf 



(9) 



where 



Ci = 2j(l - /a) - 2, 

C2 = -^(l-./a)-2, 
./a 

ML = + /a(M^,„ - Mij, 

and where Afacc, Mdon indicate binary component - donor 
and accretor masses, and the indices i, / denote the initial 
and final values, respectively. For our standard model, we 
assume that half of the mass lost from the donor is also 
lost from the system (1 — /a = 0.5) with specific angular 
momentum j = I. 

For MS and HMS donors, we follow the MT episode in 
small (< 1%) mass increments and we update the stellar 
masses and radii as well as the orbital separation and the 
Roche-lobe radii. Mass transfer is terminated when, for 
both stars, the Roche-lobe radii exceed the stellar radii. 
We also take into account possible rejuvenation of the ac- 
cretor (as described in HPT) and update its radius and 
core mass. If during this mass transfer phase the donor 
mass decreases below 0.5 Mq for MS stars or below 0.3 Mq 
for HMS star, we terminate our calculations for the spe- 
cific system, since such low-mass stars will not end their 
evolution in a Hubble time. For donors beyond the MS, 
we assume that the entire stellar envelope is lost and the 
final donor mass is equal to its core mass: He core or CO 
core mass for H-rich or He-rich stars, respectively. 

The mass-ratio criterion for stable mass transfer (see 
eq.^ also accounts for cases where the accretor is a NS 
or BH (see Kalogera & Webbink 1996; Kalogera 2000). 
We extend this treatment also to WD and allow for non- 
conservative mass transfer in cases where the criterion 
is satisfied. When WD masses exceed 1.44 M© (Chan- 
drasekhar mass according to our stellar models) we assume 
that a Type la SN occurs^ We also allow for accretion- 
induced collapse of NS to BH and we treat the boundary 
-^^max ^ model parameter. It is clear that our treat- 
ment of dynamically stable accretion onto compact objects 
does not include some important aspects (e.g., Eddington- 
limited accretion, transient systems), but it turns out that 
these MT phases are not crucial to double compact object 
formation. 

Conservative Mass Transfer. In general, we assume that 
dynamically stable mass transfer is non-conservative and 
we include the case of conservative mass transfer (no mass 
or angular-momentum loss from the system) in our param- 
eter study. In this case, /a = 1 and equation |9| reduces to 
(e.g. Verbunt & Van den Heuvel 1995), 



A,, 



M^^Ml,, 



mLM: 



f 

acc 



(10) 



*For bare CO cores, in the absence of any published models, we adopt a fixed radius of 0.01 Rq. We find that variation of this assumed 
value does not affect our results in any appreciable way (even for an increase of this radius by an order of magnitude). 

^The issue of progenitors of Type la SN is still not resolved (e.g., Branch et al. 1995). However, our assumption of the occurrence of Type 
la SN does not infiuence our results. Since accreting WD companions are less massive or of comparable mass, even if WD were to collapse to 
NS due to accretion, the chance of forming double compact objects is vanishing 
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Standard Common Envelope Phase. In all cases where 
criterion 2.7 is not satisfied, the donor has reached the gi- 
ant branch, and the accretor is either a MS, HMS, WD 
or a high mass He-rich giant branch star, we expect mass 
transfer to be dynamically unstable and a common enve- 
lope spiral in and ejection to occur (Iben & Livio 1993). 
We assume that, during the CE phase, the donor loses its 
entire envelope and the spiraling-in companion does not 
accrete any of the envelope material. 

We base our treatment on the energy formalism of 
the CE evolution (Webbink 1984), where the envelope is 
ejected on the expense of the binary orbital energy, and a 
tight post-CE system forms. The final orbital separation 
Af is calculated using. 



etc 



2Af 



2A 



don 



'don,rl 



(11) 

where, Mdon and Mace are masses of the donor and its 
companion, Mdon.cnv is mass of donor's envelope, i?don,ri 
is the Roche lobe radius of the donor, and the indices z, / 
denote the initial and final values, respectively. Parame- 
ter A describes the central concentration of the giant (de 
Kool 1990; Dewi & Tauris 2000). The right hand side of 
equation |ll| expresses the binding energy of the donor's 
envelope, the left hand side represents the difference be- 
tween the final and initial orbital energy, and ace is the CE 
efficiency with which orbital energy is used to unbind the 
stellar envelope. If the calculated final binary orbit is too 
small to accommodate the two stars then a merger occurs. 
In our calculations, we combine acc and A into one CE 
parameter, and for our standard model, we assume that 
ace X A = 1.0. 

Double Common Envelope Phase. Brown (1995) sug- 
gested that if two binary components are giants with con- 
vective envelopes and mass transfer is initiated, then a 
double CE can form, where the two stellar cores spiral in 
the combined stellar envelopes. As in the case of low-mass 
helium stars, hydrodynamical studies of such a phase have 
not been undertaken, and we note that, despite its plau- 
sibility, the validity of such an assumption has not been 
demonstrated in detail (see also footnote in §2.1.2). We 
use an analogue of equation 11 to describe the energy bal- 
ance of a double CE and to calculate the change of binary 
separation 



compact object is much smaller than the mass of the giant. 
However, for a large fraction of CE events this assumption 
is not justified, e.g., binaries with NS and a low-mass He- 
rich giant, or with a ~ 5 — 10 BH and a relatively 
massive giant companion. 

In StarTrack, we relax the simplifying assumption made 
by Bethe & Brown (1998) and derive a formulation that 
can be applied to CE phases with hyper-critical accretion 
for any compact object and companion masses. This new 
formulation involves the numerical solution of a set of or- 
dinary differential equations for the final compact object 
mass and orbital separation of the post-CE system. This 
non-approximate solution, results in tighter post-CE sys- 
tems and lower mass compact objects compared to the 
analytic expression of Bethe & Brown (1998). We present 
the detailed derivation in Appendix A. 

Unmodeled cases. There are two cases of binary interac- 
tions that are not covered by the above MT prescriptions, 
as we cannot readily foretell the outcome of some mass 
transfer events. We do not follow the subsequent evolu- 
tion of binaries going through such events, but we keep 
a record of their occurrences. For both types of unmod- 
eled MT events, the fraction of primordial binaries that 
encounter such events is below 0.1% for all our models. 

These two types of unmodeled events are: (i) Unevolved 
donor with an evolved companion (i.e., giant with con- 
vective envelope). This configuration is encountered very 
rarely, as usually it is the giants that fill their Roche lobes, 
(ii) Compact objects accreting from massive evolved he- 
lium stars (with radiative envelopes). We expect that most 
of these cases lead to stellar mergers. Moreover, if any two 
of MS, HMS or massive HGB stars fill their Roche lobes 
at the same time, the formation of a contact system is 
recorded, and the calculations for these systems are termi- 
natedfl. 



GM(mI GMIM^\ _ GMlMi^cnv GM|Af2,c 



2At 



2Ai 



XR 



2,rl 



where Mi,M2 are the two donor masses. Mi, 



,Mo 



(12) 



2.2.5. Core Collapse and Supernova Explosion Events 

For each massive star, the time of core collapse is set by 
the single star models (taking though into account mass 
variations due to winds and binary interactions). When 
either component of a binary reaches this stage, we gen- 
erate a random location in the orbit for the event to take 
place (note that for eccentric binaries this choice will af- 
fect the outcome, since the components' separation and 
relative velocities are different at different locations in the 
orbit). The core-collapse event is assumed to be instan- 
taneous and the mass of the remnant is calculated using 
equation |l|. Note that if the remnant is formed through 
complete fall back (leading always to direct BH formation) , 
we do not expect a SN explosion (hence no kick and no 
mass loss) and the orbit remains unchanged (Fryer 1999). 
When BH is formed through partial fallback we treat the 
event as a SN explosion (see Podsiadlowski et al. 2001). 

We calculate the effect of a SN event on binaries in three 
steps: (i) We estimate the mass of the remnant. The rest 
of exploding star is immediately lost from the binary (with 
the angular momentum specific to the exploding compo- 
nent). We assume that the ejecta do not have any effect 
on the companion (e.g., Kalogera 1996). (ii) We calculate 

®We find that the fraction of contact systems formed is rather independent of the adopted evolutionary model and is in the range 4—6% of 
the total simulated binary population. 



are their envelope masses and other symbols have the same 
meaning as in equation pT| (see also Nelemans et al. 2001a). 

Common Envelope with Hyper- critical Accretion. It has 
been suggested that if a NS or BH evolves through a CE 
phase, the compact object may accrete significant amounts 
of material because of hyper-critical accretion (Blondin 
1986; Chevalier 1989, 1993; Brown 1995). Bethe & Brown 
(1998) derived an analytic scheme for the evolution of the 
binary orbit and the increasing mass of the accreting com- 
pact object, under the assumption that the mass of the 
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the compact object velocity which is the vector sum of the 
orbital velocity of the pre-coUapse star at the orbital posi- 
tion and the kick velocity. The kick velocity is assumed to 
be randomly oriented and its magnitude is drawn from an 
assumed distribution. The kick magnitude is also scaled 
with the amount of material ejected in the SN explosion, 

^kick - (1 - fib) * V, (13) 

where V is the kick magnitude drawn from a given as- 
sumed distribution, and /fb is a fall back parameter de- 
fined in §2.1.3, and the Vkick is the kick magnitude we 
use in our calculations. For NS remnants and no fall back 
(/fb — 0) and Vkick = V". In our standard model we use 
a kick magnitude distribution very similar to the one de- 
rived by Cordes & Chernoff (1998): a weighted sum of 
two Maxwellians, one with a = 175 km s~^ (80%) and 
the second with a = 700 km s"^ (20%). This distribu- 
tion accounts for the fact, that besides average velocity 
pulsars a significant fraction of neutron stars have veloci- 
ties above 500 kms~^ (Arzoumanian, Cordes, & Chernoff 
1997). (iii) We calculate the total energy (potential and 
kinetic) of the new orbit for the remnant (new velocity 
and mass, same relative position) and its companion. If 
the total energy is positive, then the system is disrupted, 
and its components will evolve separately. We calculate 
their subsequent evolution as single stars, but we do not 
follow their trajectories in the present study. If the total 
binary energy is negative, the system after SN explosion is 
bound, and we calculate its new parameters (e and A). We 
also check whether the two components have merged due 
to the SN mass loss and kick in which case we terminate 
the evolution. Finally, we calculate the post-SN center of 
mass velocity of a binary. 

2.2.6. Orbital Decay Due to Gravitational Radiation 

Once both components of a binary have ended their nu- 
clear evolution and formed stellar remnants, their orbits 
change only due to angular momentum loss through gravi- 
tational radiation. We calculate the merger times of these 
binaries using the formalism developed by Peters (1964), 
the result of which is briefly described below. 

The decay rates of eccentricity e and orbital separation 
A are combined to derive the eccentricity evolution with 
time e{t) and integrate to the limit e — > (as the system 
approaches merging, i.e., ^ ^ 0): 

en 

_ 12 f e^9/i9[i + (I21/304)e^]"si/^^99^e 



where cq and Aq are the initial eccentricity and orbital 
separation. 



" eoi2/i9[H- (121/304)eo2]870/2299' 

64 G^MilVhiMi + Ah) 

and where G is the gravitation constant and c is the speed 
of light. In practice the limit A — > (accompanied by 



e — > 0) is unrealistic, since a merger occurs earlier, but 
the very final parts of this integration do not contribute 
significantly to the merger time, as the very late stages of 
inspiral proceed on very short timescales. For example, 
the merger time we calculate for PSR B1913+16 is within 
1% of the result obtained by solving the system of two or- 
dinary equations for the eccentricity and orbital separation 
(e.g.. Junker & Schaefer 1992). 

We note that, for a circular orbit, we obtain a merger 
time: 

W-g = 4/{Ap). (15) 

2.2.7. Star Formation History and Binary Fraction 

For most of our models, we assume that star forma- 
tion has been continuous in the disk of our Galaxy for the 
last lOGyr as inferred from observations (e.g., Gilmore 
2001) and as predicted by recent theoretical modeling (e.g., 
Kauffmann, Chariot & Balogh 2001). We start the evolu- 
tion of a single or a binary system tbirth ago, and follow 
it to the present. The birth time tbirth is drawn randomly 
within the range 0-10 Gyr. In our parameter study, we 
examine a model of instantaneous star formation. 

For most of our models, we also assume a binary frac- 
tion of /bi = 0.5, which means that for any 150 stars we 
evolve, we have 50 binary systems and 50 single stars. We 
treat /bi as a model parameter and include it in our pa- 
rameter study. Finally, in all our models we assume solar 
metallicity Z = 0.02. 

3. RESULTS 
3.1. Population Synthesis Parameter Study 

We perform an extensive parameter study in order to 
assess the robustness of our population synthesis results. 
In Table |l| we summarize all our models based on those 
assumptions that are different from our standard model 
(model A). Parameter values and adopted distributions in 
the standard model are as given in the two previous sec- 
tions: §2.1 and §2.2. 

In models marked with letter B, we use different distri- 
butions of SN kick magnitudes imparted to NS at birth. 
In model BI we assume symmetric SN explosions, whereas 
in models B2-12 we draw kick velocities Vk from a single 
Maxwellian: 

5(14) cxl4'exp[~(Vfe/(7)2] , (16) 

varying a values in the range 10 — 600kms~^. In model 
B13 we use a kick distribution of the form suggested by 
Paczynski (1990): 

/(t4)«[l + (yfe/a)2]-i, (17) 

which allows for a significant fraction of low-magnitude 
kicks. We use a = 600km s~^, which gives a reasonable 
fit to the population of single pulsars in the solar vicinity 
(Hartman 1997). 

In model C, compact objects (NS or BH) are not allowed 
to accrete any material in CE events. In models Dl-2, we 
reduce our conservative maximum limit on the NS mass of 
3 Mq down to 2 and 1.5 Mq. Models El-3 present evo- 
lution with different effective CE efficiencies (ace x A). In 
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model Fl, we significantly decrease the amount of material 
accreted by companions in non-conservative MT events, 
whereas in model F2, we examine the case of conservative 
mass transfer. In models Gl-2, we vary the wind mass loss 
rates. In model Gl, we decrease the wind mass loss rate 
by a factor of 2, for all stars and at all evolutionary stages, 
whereas in model G2, we enhance it by factor 2. In other 
words, we calculate the wind mass loss rate as described 
in § 2.1.1, and multiply it by /wind given in Table |l|. Model 
HI corresponds to evolution with Afconv = 4.0Mq, the up- 
per mass limit for helium giants with convective envelopes 
that would initiate a CE phase. Model 112 does not allow 
for CE evolution for any Roche-lobe filling helium stars 
regardless of their mass (effectively Mconv = 0). In model 
I, the star formation rate (SFR) is altered, and instead of 
continuous we assume a burst-like star formation history, 
with all stars being formed 10 Gyr ago. In model J, we use 
an initial mass function (IMF) with a Salpeter (1955) ex- 
ponent of —2.35, and in models Kl-2, we vary the binary 
fraction of binary systems in the primordial stellar popu- 
lations. In models LI and L2, we vary the specific angular 
momentum of material lost during dynamically-stable MT 
events. In models Ml and M2, we use different distribu- 
tions of initial binary mass ratios. Model N corresponds to 
the artificial case of evolution without allowing for the ex- 
pansion of evolved helium stars. This model is nonphysical 
since low mass helium stars are known to evolve to giant- 
like stages and to significantly increase their radii (for a 
detailed discussion see Belczynski & Kalogera 2001 and 
references therein). We include this model in our param- 
eter study only for the purposes of comparing our results 
to earlier studies that did not not account for helium star 
radial evolution. In model O, we increase the mass range 
of compact object formation through partial fall back, re- 
sulting to no BH forming through a direct collapse. 

3.2. Statistical Accuracy 

Apart from the systematic uncertainties of our results 
which we examine with the broad parameter study dis- 
cussed above, we also examine the statistical accuracy of 
our population synthesis models, to make sure that we do 
not burden the quantitative predictions with unnecessary 
uncertainties. To determine the intrinsic statistical ac- 
curacy of our synthesis models, we performed 30 different 
realizations of our standard model (A), each with N — 10^ 
independently generated primordial binaries. We decide to 
limit ourselves to iV = 10^ runs, given our computational 
resources. 

In Table |2| we present the results we obtained by com- 
paring the 30 different runs. For each of the three DCO 
populations, we list the mean number of systems formed 
from the set of iV = 10^ runs, its standard deviation (as a 
percentage), as well as the maximum variation in the pre- 
dicted numbers of systems. We find that the standard de- 
viations do not exceed 10% for any type of DCO, whereas 
the maximum variation remains below 10% for the most 
frequent NS-NS binaries and below 30% for the least fre- 
quent population of BH-NS binaries. In what follows, we 
present results for our standard model using the combined 
total of these 30 runs (or their subgroups) and therefore 
statistical inaccuracies are greatly reduced and become in- 
significant compared to the systematic uncertainties. For 



some of the other models in our parameter study (e.g.. El, 
G2) we also considered results from runs with more than 
10^ primordial binaries to ensure that our statistical errors 
remain well below 10% for all three DCO types. 

3.3. Formation Paths of Double Compact Objects 

We consider double compact objects with NS or BH 
(NS-NS, BH-NS or BH-BH binaries) with merger times 
shorter than 10^*^ yr, thus coalescing DCO. In our standard 
evolutionary model, the population of coalescing DCO is 
dominated by NS-NS systems (61%), with a significant 
contribution by BH-BH binaries (30%), and a small contri- 
bution by BH-NS objects (9%). In what follows we discuss 
the main qualitative characteristics of the important for- 
mation paths (with relative formation frequencies higher 
than 1%) as well as their origin. 

In Table | we present the most important formation 
channels of coalescing DCO, for our standard model. For- 
mation channels of NS-NS, BH-NS and BH-BH binaries 
are marked by NSNS, BHNS and BHBH, respectively, 
they are listed in order of decreasing relative formation 
frequency (second column) with respect to the whole DCO 
coalescing population. The details of each evolutionary se- 
quence, i.e., MT episodes and SN explosions are also given. 
Results were obtained based on the evolution of 3 x lO'' 
primordial binaries. 

3.3.1. Populations of Double Neutron Stars 

Using the StarTrack population synthesis code, we iden- 
tified a number of new NS-NS formation channels (see be- 
low) . This is a result of two improvements in the imple- 
mentation of our population synthesis code, since our pre- 
vious short study of the NS-NS formation (Belczynski & 
Kalogera 2001, hereafter BKOl). First, we have replaced 
the approximate prescription suggested by Bethe & Brown 
(1998) for the hyper-critical accretion during CE phases, 
with a newly derived numerical solution (see Appendix). 
Second, we allow for hyper-critical CE evolution of low- 
mass helium stars with compact objects. In the studies 
presented in BKOl, we allowed binaries with low-mass he- 
lium giants to evolve through DCE and standard SCE, 
but we had assumed that CE events of helium giants with 
compact objects lead to mergers, and possibly a gamma- 
ray burst (e.g.. Fryer et al. 1999). However, due to the 
small mass of helium giant envelope at the onset of CE 
event (^ 1 — 1.5Mo), we find that these systems survive 
the CE events, and form very tight NS-NS binaries (see 
also footnote in §2.1.2). 

Double neutron stars are formed in various ways through 
more than 14 different evolutionary channels identified in 
Table |[ However, the entire population of coalescing NS- 
NS systems, may be divided into three subgroups. 

Group I. This subpopulation consists of non-recycled 
NS-NS systems, first identified by BKOl. These are sys- 
tems in which none of the two NS ever had a chance of 
getting recycled through accretion. Our current results for 
the predicted formation rates and properties of the non- 
recycled NS-NS systems, have not been affected by the 
two improvements discussed above. As shown in Table ^ 
these recently identified non-recycled NS-NS systems are 
formed via the NSNS:09, NSNS:11 and NSNS:12 channels. 
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NS-NS systems (belonging neither to Group I nor II) 
formed, through more or less classical channels (Bhat- 
tacharya & van den Heuvel 1991). The formation path 
denoted NSNS:14 corresponds to what is usually consid- 
ered to be the "standard" NS-NS formation channel (Bhat- 
tacharya & van den Heuvel 1991). Since we account for 
hyper-critical accretion in CE, the formation rate is de- 
creased because some NS (but not all, as assumed by 
Portegies Zwart & Yungelson 1998 and by Fryer et al. 
1999) collapse to BH. Furthermore our treatment of the 
hyper-critical accretion typically leads to tighter post-CE 
systems, causing more binaries to merge in CE events, and 
thus decreases the number of possible NS-NS progenitors. 

Group II strongly dominates the population of coalesc- 
ing NS-NS systems (81%, for standard model calculation) 
over group III (11%) and I (8%). In general this character- 
istic is preserved in all, except a few extreme evolutionary 
models. The discussion of the dependence of the formation 
rates of the various NS-NS subpopulations on population 
synthesis model assumptions will be presented in Belczyn- 
ski et al. 2002a. 

3.3.2. Populations of Black Hole Binaries 

In general, BH-NS and BH-BH binaries are formed 
through just a few distinct channels, with a moderate num- 
ber of MT events (2-3), in contrast to our findings for 
NS-NS systems. 

Helium star evolution, radial expansion and CE phases 
are much less important for the formation of BH-NS and 
BH-BH binaries. The reason is that for most of these pro- 
genitors the first-born compact object is massive enough 
that even when helium stars evolve to the giant branch, 
they do not expand to large radii nor they lead to pos- 
sible CE evolution (see MT criteria in § 2.2.4 and chan- 
nel BHNS:03 in Table |). Instead these DCO form most 
efficiently through channels that closely resemble those 
NS-NS conventionally thought to be "standard" (Bhat- 
tacharya & van den Heuvel 1991; Fryer et al. 1999): evo- 
lution is initiated with a phase of non-conservative mass 
transfer and followed either by a CE phase or the forma- 
tion of the first compact object (see BHNS:01, BHNS:02, 
and BHBH:01, BHBH:02). 

3.4. Physical Properties of Double Compact Objects 

3.4.1. Double Compact Objects Component Masses 

Standard Model. In Figure E 
tributions of DCO components 

with first- and second-born compact objects in the top and 
bottom panel, respectively. 

The shape of the calculated mass distributions can be 
understood as a convolution of the initial-remnant mass 
relation for single stars (top panel of Figure □) with the 
adopted initial stellar mass function (see eq. g). Devia- 
tions from such a convolution reflect the effects of binary 
evolution on the compact object masses. 

First, we discuss the mass distribution of the second- 
born compact objects, as they are less affected by binary 
evolution. This distribution starts at the lowest possible 
NS mass (allowed by our single star models) of 1.2M0, 

'^Mass distributions of compact objects with white dwarf companions as well as single compact objects (formed both from single star 
progenitors and from components of disrupted binaries) are presented in Belczynski, Bulik &; Kluzniak 2002b. 



which involve DCE of two low-mass helium giants, which 
were already allowed in the earlier version of StarTrack. 

The unique qualitative characteristic of this NS-NS for- 
mation path is that both NS have avoided recycling. The 
NS progenitors have lost both their hydrogen and helium 
envelopes prior to the two supernovae, so no accretion from 
winds or Roche-lobe overflow is possible after NS forma- 
tion. Consequently, these systems are detectable as ra- 
dio pulsars only for a time (~ 10® yr) much shorter than 
recycled NS-NS pulsar hfetimes (~ 10^ - 10^° yr in the 
observed sample). Such short lifetimes are of course con- 
sistent with the number of NS-NS binaries detected so far 
and the absence of any non-recycled pulsars among them. 

We note that the identification of the formation path for 
non-recycled NS-NS binaries stems entirely from account- 
ing for the evolution of helium stars and for the possibility 
of double CE phases, both of which have typically been ig- 
nored in previous calculations (with the exception of Fryer 
et al. 1999, where, however, such events were assumed to 
lead to mergers). 

Group II. This subpopulation consists of tight, short 
lived binaries with one recycled pulsar. Their merger times 
are typically ~ 1 Myr or even smaller (see §3.4.5). As 
shown in Table ||, these new dominant NS-NS systems are 
formed via the NSNS:01-08, NSNS:10 and NSNS:13 chan- 
nels, with the common characteristic that the last binary 
interaction is a hyper-critical CE of a low-mass helium gi- 
ant and the first-born NS. 

In Belczynski, Bulik, & Kalogera 2002a we describe in 
detail the formation of a typical NS-NS binary of group II. 
The most channel identified as the most efficient for NS- 
NS formation (NSNSiOl) corresponds to the "standard" 
channel of Fryer et al. (1999). The only difference is an 
extra CE event which originates from allowing for helium 
star evolution and without a priori assumptions about the 
CE outcome. The second most dominant channel, involv- 
ing two consecutive MT episodes and then two SN explo- 
sions, closely resembles our channels: NSNS:02, NSNS:04, 
NSNS:06, NSNS:10, NSNS:11, NSNS:12. The only differ- 
ence again remains an extra MT episode from evolved, 
Roche-lobe-filling helium stars. 

The most dramatic effect of the binary evolution up- 
dates is reflected in the existence of a whole new popu- 
lation of coalescing NS-NS stars formed in the Group II. 
In our standard model these channels contribute 50% of 
the DCO population, and their common characteristic is 
that the last binary interaction is a hyper-critical CE of a 
low- mass helium giant and the first-born NS. It turns out 
that the majority of these systems survive the HCE event 
and form tight NS-NS binaries. Had we not taken into ac- 
count the radial expansion of low-mass helium- rich giants, 
the progenitors of this dominant NS-NS population would 
have evolved without any further MT. Most of them would 
have still formed NS-NS systems, however, not as tight as 
after this last CE episode. We have actually examined this 
alternative and found that about half of them would have 
formed binaries with merger lifetimes longer than 10^" yr. 
Once again, we see the importance of helium star evolution 
on DCO population synthesis. 

Group III. This subpopulation consists of all the other 



we present the mass dis- 
for our standard model, 



11 



rises sharply with a peak at « 1.4M0 and then dechnes 
down to w 3M0 (the assumed M^f^ for our standard 
model). This strong peak at low masses is the result of 
a rather weak dependence of NS mass on ZAMS mass 
(Figure |l|) combined with the steep IMF. Following the 
initial peak, the distribution becomes flat in the range 
« 3 — IOM0. This flattening is caused by the balanced ef- 
fects of a slightly rising initial-remnant mass relation and 
a declining IMF. 

At ~ IO.SMq the distribution peaks again, and then 
rapidly declines. This final peak corresponds to the sat- 
uration of BH masses at IO.SMq for a wide range of 
ZAMS masses in the initial-remnant mass relation. 

For the model described here, the maximum mass of 
the second-born compact objects is ~ ISM©, in excess of 
the maximum mass of w IIMq of compact objects formed 
in single star evolution (Figure [l]). Such slightly more 
massive compact objects are allowed only because of MT 
episodes involving their progenitors. In some cases they 
can accrete enough material to significantly increase the 
core mass (rejuvenation through accretion), and form more 
massive remnants. 

The mass distribution of first-born compact objects 
shows more prominently the effects of binary evolution, 
since the compact objects themselves and not just their 
progenitors might have been affected by mass transfer 
events. The three basic qualitative features are still 
present, with an initial peak, a subsequent flattening, and 
a smaller final peak. However, there are two easily iden- 
tified differences between the two distributions. First, the 
initial peak is lower and broader, and second, the peak at 
high masses shifts to « 13Mq, and the maximum compact 
object mass increases to w 14M0. Both effects are due to 
accretion onto the first-born compact objects from their 
non-degenerate companions. 

Our standard-model mass distribution of compact ob- 
jects in coalescing double compact objects peaks at 
1.4M0 for NS and then extends up to 14M0 for BH. 
This agrees well with the observed NS masses in binary 
stars: 1.1 - I.6M0 (e.g. Thorsett & Chakrabarty 1999) 
and also with estimated masses of black hole binary can- 
didates: ~ 3 — 2OM0 (for references see §2.1.3). It is im- 
portant to note that the maximum BH mass of ~ 14M0 
can be increased up to ~ 20 — 25M0. This can be achieved 
either by decreasing the wind mass loss rate by factor of 2 
(see the bottom panel of Figure |l|) and/or by increasing the 
amount of material accreted by stars in non-conservative 
MT events. 

We can also compare results of our population synthesis 
calculations to the theoretical compact object mass dis- 
tributions presented by Fryer & Kalogera (2001, hereafter 
FKOl; their Figures 7 and 8). We find some striking sim- 
ilarities. The single star initial- remnant mass relations 
are very similar as expected given that our calculation of 
compact object masses are based on the same core-collapse 
hydrodynamical calculations (Fryer 1999). FKOl find that 
in the range of compact object formation 80% of their sin- 
gle star remnants are NS, and the rest are BH. We find 
that for the standard model of single star evolution 81% 
of remnants are NS, and 19% are BH. 

More importantly, the shape of the FKOl CO mass dis- 
tribution formed in the presence of stellar winds and bi- 



nary companions, resembles our distribution, even though 
FKOl did not use population synthesis. Apart from the 
qualitative similarities, we also note that our results con- 
firm the findings of FKOl for a continuous mass distribu- 
tions and the absence of a gap or isolated narrow spike at 
about 7M0, claimed by Bailyn et al. (1998). Recent ob- 
servational results also point in the direction of a broader 
range of BH masses (Froning & Robinson 2001; Orosz et al. 
2001). However, there are still two significant differences 
between our results and those of FKOl. One is related to 
the maximum BH mass (about lOM© in FKOl compared 
to our I4M0) and the other to the contribution of BH to 
the compact object population. Both are linked to a num- 
ber of binary interactions that are not taken into account 
by FKOl and have the effect of increasing the masses of 
compact objects. In our standard-model calculations we 
find that, in double compact objects, NS and BH repre- 
sent 65% and 35% of the compact objects, respectively, 
because of accretion effects. In contrast, for the case in 
FKOl where effects of stellar winds and binary compan- 
ions, the contribution of BH actually decreases compared 
to single star calculations. 

Parameter Study. We have explored a large number of 
population models and have found that the main qualita- 
tive features of the mass distributions remain unaffected: 
the low-mass peak from NS formation and the flat form of 
the distributions for a wide range of BH masses are present 
in many different models (e.g., Bl, B12, F2, J). Accretion 
effects and the less prominent peak at higher masses are 
also clearly seen in almost all of the models. 

Nevertheless some dramatic quantitative differences are 
revealed by a few models. In the case that all wind mass- 
loss rates are increased by a factor of 2 (G2), the maximum 
compact object mass turns out to be only ~ 3.5 Mq and 
therefore such a case appears to be highly unlikely. On 
the other hand, the model where all wind mass-loss rates 
are decreased by a factor of 2 (Gl) allows for BH masses 
as high as ~ 21 Mq. 

Another class of models that show significant quantita- 
tive differences from the standard model are those with 
varying CE efficiencies, especially those with low values. 
The results for model El (acE x A = 0.1) are shown in 
Figure ^j. Although the basic shape of the distributions 
resembles that of the standard model, the second-born 
compact objects span a much narrower range in mass: 
1.15 — 9.7 M0 (with most of compact objects formed with 
mass smaller than GMq). Low CE efficiency discrimi- 
nate against tight pre-CE systems, cause a large number 
of DCO progenitors to merge, and therefore DCO rates 
to decrease significantly. In particular, the effects appear 
to be most prominent for BH-BH binaries where massive 
second-born BH are essentially eliminated. As shown in 
§ 3.3 and Table ||most BH-BH progenitors experience first 
a non-conservative MT and then a CE episode. For the 
case of low CE efficiencies, CE survival is favored for rather 
wide pre-CE systems with donors (stars that initiate the 
CE episode and progenitors of the second-born compact 
objects) with low envelope masses. In these cases CE ejec- 
tion is facilitated, but at the same time lower mass pro- 
genitors are favored, leading to a restricted mass range for 
these second-born compact objects. At the same time the 
bias in favor of wide pre-CE systems translates to a bias in 
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favor of progenitors where the first MT episode is initiated 
late in the evolution of the primary, when the stellar radii 
are larger and most of the stellar envelope has burnt into 
a high mass core. This bias clearly favors the formation of 
higher mass first-born compact objects, leading overall to 
BH-BH binaries with rather extreme mass ratios. 

Given the sensitivity of the DCO mass distribution to 
the CE efficiencies, it may be possible in the future to dis- 
criminate between models based on mass measurements 
of BH-BH binaries from gravitational-wave detections and 
data analysis. 

3.4.2. Double Compact Objects Orbital Separations 

Standard Model. The distribution of DCO over orbital 
separations at the time of their formation is shown in Fig- 
ure ^ for our standard model. 

The distribution for the whole population (top panel) 
covers a wide range of values ~ 0.1 — IOOR0 and is char- 
acterized by two distinct peaks. It is evident from the 
bottom panel that the peak at shorter separations is due 
to NS-NS binaries and the one at wider is due to BH bi- 
naries. Examination of the NS-NS distribution indicates 
that the majority of NS-NS systems, are formed with very 
small orbital separations (< IR©), but there is a long tail 
stretching to much larger orbital separations (even beyond 
IORq). This result is related to our discussion of the main 
NS-NS formation channels (§3.3). The NS-NS popula- 
tion consists of three distinctive groups (see § 3.3.1), which 
shape the orbital separation distribution. "Classical" NS- 
NS systems (group HI) that contain recycled pulsars have 
wide separations, mostly in the range ~ 1 — lOR© (where 
the two observed systems lie). However, NS-NS binaries in 
which none of the NS ever had a chance of getting recycled 
(group I) and very tight NS-NS binaries formed via the 
newly identified channels (group II) actually dominate the 
produced NS-NS population and produce a strong peak at 
much lower orbital separations. These differences in or- 
bital separations affect the derived merger times and have 
important implications for the detection of such systems, 
which we discuss in §3.5.2. 

Typical BH binaries have much wider separations than 
NS-NS systems (Figure E|) . This is a result of the lower (or 
zero) SN kicks impartea to BH which allows the survival 
of wider progenitors, typical of the more massive binaries 
required for the formation of the BH binaries. Note that 
BH-NS binaries tend to have intermediate size orbits, af- 
fected in part by the high kicks imparted to one of the two 
compact objects. 

Parameter Study. The double-peaked shape of the or- 
bital separation distributions found in the standard model 
persists in a robust way in the majority of the models in 
our parameter study. The full range of separation values 
also appears to be robust (0.1-100 R©). However, the rela- 
tive height of the two peaks varies significantly from model 
to model, primarily because of the varying contribution of 
different types of systems to the DCO population. The 
reasons and a discussion of these variations are given in 
§3.5.1. 

Here, we identify a couple of extreme cases: In the mod- 
els with zero kick velocities (Bl) or increased wind mass- 
loss rates (G2), BH binaries have a very small formation 
rate relative to NS-NS systems (see also Portegies Zwart & 



Yungelson 1998; Fryer et al. 1999), and therefore the NS- 
NS peak at shorter orbital separations strongly dominates. 
In contrast, models with low CE efficiencies hamper the 
formation of DCO in tight orbits and the peak due to BH 
binaries strongly dominates. In model H2, as expected, the 
peak at low separations for NS-NS disappears altogether, 
since no helium-star CE episodes are allowed. 



3.4.3. Double Compact Objects Eccentricities 

Standard Model. In Figure |^, we present eccentricity 
distributions of coalescing compact object binaries for our 
standard model. The values at the time of DCO formation 
are quite high, e > 0.1 for the majority of the population, 
and of course originate from the fact that the last stage 
prior to DCO formation typically involves a SN explosion 
and asymmetric kicks. The fraction of DCO formed in cir- 
cular orbits is rather small, and occurs when the second 
compact object is formed through a direct collapse to a BH 
(BH-BH: ~ 2.4% and BH-NS:~ 0.1%). Portegies Zwart & 
Yungelson (1998) and Fryer et al. (1999) have also found 
that NS-NS population is born with high eccentricities, 
which is in good agreement with our results. 

Apart from the eccentricities at formation, we also ex- 
amine the distribution at later stages when gravitational 
radiation has caused the binary orbit to shrink and the 
gravitational- wave frequency has increased to ~ 40 Hz, 
the lower end of the LIGO I band. So far DCO have 
been assumed to be circular for purposes of gravitational- 
wave detection and data analysis. We examine whether 
our model population satisfies this assumption due to an- 
gular momentum losses and circularization. We use the 
dependence of orbital separation on eccentricity due to 
gravitational-wave emission as derived by Peters (1964): 
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where cq is given in equation ^ and depends on the orbital 
separations and eccentricities at DCO formation. We use 
our results for the masses to calculate first the orbital sep- 
aration at 40 Hz and the corresponding frequency. The 
resulting distribution is shown in Figure ^ (solid line). 
It is clear that by the time these coalescing DCO en- 
ter the LIGO I band, the eccentricities are small enough 
(e < 0.0001) that the assumption of circular orbits is well 
justified. 

It is interesting to note that, although the eccentricity 
distribution at birth is single-peaked, the distribution at 
later stages becomes double-peaked. This is a direct re- 
sult of the double-peaked form of the orbital separation 
distribution. In more detail, the peak at lowest eccentrici- 
ties is populated by the heavier BH binaries and the peak 
at somewhat higher eccentricities by NS-NS binaries. The 
distinction is related to the weaker effects of gravitational- 
wave emission for the less massive systems. 

Parameter Study. The shape and range of typical val- 
ues for the orbital eccentricities appear to be quite robust. 
Differences becomes noticeable only in models of varying 
kick-magnitude distributions. As expected, the typical ec- 
centricities at formation decrease with a decreasing aver- 
age kick magnitude. 
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The variations in the shape and position of the eccentric- 
ity peaks for later stages are more prominent and closely 
follow the corresponding variations in the distributions 
over orbital separations. Nevertheless, the eccentricities 
at 40 Hz remain below 10^'^ for all the models we have 
examined. 



3.4.4. Double Compact Objects Center of Mass Velocities 

Standard Model. Binaries are expected to acquire sys- 
temic (center-of-mass) velocities after core-collapse events 
because of the combined effects of mass loss and SN kicks 
imparted to the compact remnants. Understanding of 
these velocities is crucial in studies of dynamical evolu- 
tion of these populations, and their kinematic properties 
and spatial distributions with respect to host galaxies. The 
distributions of center-of-mass velocities after the first and 
second core-collapse events are shown in Figures | and 0, 
respectively. Top panels on both figures show the distri- 
bution for the whole DCO population, and bottom pan- 
els show the distributions separately for each of the DCO 
groups. 

It is evident that there are major differences between 
the two plots. The distribution after the first core-collapse 
event (Vi) shows a well defined narrow peak at low veloc- 
ities (~ 25km s" ) and a slowly decaying tail at higher 
velocities. After the formation of the second compact ob- 
ject velocities (V2) populate a much broader peak posi- 
tioned at high magnitudes (~ 200km s~^) with a velocity 
tail reaching to very high velocities (> 600kms~ ). More- 
over, in these distributions a narrow spike at 0km s~ is 
also identified. 

All these characteristics can be naturally explained by 
the different binary orbital properties of binaries at the 
time of the first and the second compact object forma- 
tion. For the majority of DCO progenitors, the first core- 
collapse event happens when the binary orbit is still wide, 
accommodating two massive stars and is affected primar- 
ily by mass loss through winds and non-conservative mass 
transfer events. The typical orbital velocities in these pre- 
collapse systems are ~ 10 — 50kms~^ and it is mostly 
these orbital velocities that determine the typical post- 
collapse systemic velocities of bound systems (Kalogera 
1996). For Maxwellian kick magnitude distributions an 
upper and lower limit on the systemic velocity can be de- 
rived analytically and it is found to be independent of 
the average kick magnitude (depends only on the stellar 
masses and the pre-collapse relative orbital velocity; see 
Kalogera 1996). As a secondary effect, within the range of 
values defined by these two limits, systems tend to acquire 
higher systemic velocities when higher kicks are imparted. 
A smaller fraction of DCO progenitors experience a CE 
phase prior to the first core-collapse event, and have much 
tighter orbits (orbital velocities of order ~ lOOkms"^), 
leading to higher systemic velocities. The more subtle in- 
fluence of the average kick magnitude is evident when look- 
ing at the distributions for each DCO class, where NS-NS 
progenitors are slightly shifted to higher and BH binaries 
to lower systemic velocities (NS-NS progenitors tend to be 
lighter than BH-BH progenitors). 

At the time of the second SN explosion, the population 
of compact object binary progenitors is dominated by tight 



binaries with relative orbital velocities in the low hundreds. 
The result is a strong, broad peak at 100 — 300 km s^^ with 
a tail to higher velocities (Figure Systems with lower 
systemic velocities do form, but the majority of them are 
so wide that their merger times exceed the Hubble time. 
Another noticeable difference from Figure ^ is that the 
BH-NS population now closely follows the NS-NS popula- 
tion. This is understood as both the pre-collapse orbital 
characteristics and the typical kick magnitudes are similar 
for these two populations at the second compact object 
formation, since in the majority of BH-NS systems the NS 
is formed second. It is evident that typically BH-BH bina- 
ries acquire the lowest systemic velocities. This is a result 
of three combined effects: they are formed from wider pro- 
genitors (§3.4.2), they are heavier systems, and some BH 
are formed via direct collapse and hence do not acquire 
any kick. Strikingly, in both Figures ^ and there is a 
prominent spike at zero systemic velocity, which is popu- 
lated by BH-BH systems where both BH were formed via 
direct, assumed to be symmetric collapse. We note that, 
although these systems tend to be wider than NS-NS, they 
still merge within a Hubble time because of the higher 
masses, hence stronger gravitational radiation involved. 

It is useful to note one other implication of the small sys- 
temic velocities after the first core-collapse event. These 
combined with the fact that the time between the two col- 
lapse events is much shorter than typical coalescence times, 
implies that it is the systemic velocities after the second 
collapse (along with the merger times) that determine the 
spatial distributions of DCO merger sites with respect to 
host galaxies. 

Portegies Zwart & Yungelson (1998) presented results 
on NS-NS systemic velocities, and Fryer et al. (1999) cal- 
culated systemic velocities for both NS-NS and BH-NS 
binaries. Comparison with our results shows overall good 
agreement; both Portegies Zwart & Yungelson (1998) and 
Fryer et al. (1999) found that most systems acquire sys- 
temic velocities of the order of ^ 200 — 300 kms^^, with 
some binaries being accelerated to much higher velocities 
(- 500 - 1000 kms-i). 

Parameter Study. The distribution of systemic veloci- 
ties after the first core-collapse event is found to be very 
robust in our parameter study. The most clear, but not 
big, changes in the distribution shape are visible in models 
which affect the pre-SN orbital separations, e.g. in mod- 
els El and F2. For these models, the distribution starts 
with an initial peak, as for the standard model, however 
the high velocity tail is depleted at velocities below 100- 
150 kms~^. All of DCO progenitors experience either CE 
phase or MT event prior to the first SN explosion. De- 
creasing the CE efficiency (model El) or increasing the 
fraction of mass accreted by companions in MT events 
(model F2) lead to a decrease in post-CE or post-MT sep- 
arations. Therefore, at the time of the first SN explosion, 
progenitors have higher orbital velocities, and DCO pro- 
genitors in models El and F2 end up with higher systemic 
velocities than in the standard model, depleting the dis- 
tribution of low velocity values. Some small changes are 
also seen with kick magnitude variations. For very small 
kicks, the high velocity tail of the standard model distribu- 
tion disappears, whereas, for very high kicks, many more 
systems are disrupted, but most of the ones that survive 
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the collapse, end up with somewhat higher systemic veloc- 
ities. As a result, the low velocity peak of the distribution 
is depleted and the tail extends to higher velocities. 

Similarly, the qualitative features of the distributions af- 
ter the second compact object formation are robust, with 
an exception for the spike at zero velocity, the presence 
of which depends on the contribution of BH-BH binaries 
to the DCO population (e.g., models Bl, B6, B13, D2, 
F2, J, and the extreme case of G2 where essentially no 
BH-BH binaries are formed). For a few models, we find 
some significant differences, but are all explained by the 
relative contribution of the various DCO groups. For ex- 
ample, for very high kicks (model B12) and in the absence 
of helium-star CE phases (model H2), the formation rate of 
coalescing NS-NS is so small that the velocity distribution 
is dominated by BH-BH binaries, and in particular those 
which do not receive a kick (i.e., the zero- velocity spike 
becomes most dominant). The few surviving NS binaries 
populate a much broader range of systemic velocities with 
low normalization and a roughly fiat tail out to hundreds 
of kms~^. 



these two earlier studies. 

Parameter Study. Unlike all other binary properties 
we have examined, the qualitative characteristics of the 
merger-time distributions appear to be rather sensitive to 
a number of model parameters. This is understood in 
terms of the strong dependence of merger times on initial 
separations and eccentricities, and, although their distri- 
butions change only slightly with model parameters, the 
corresponding change of merger time distributions is quite 
dramatic. In Figure 



we show three of the most different 
DCO merger time distributions from models with some- 
what extreme assumptions. 

In the top panel of Figure H (model with zero kicks, Bl), 
the DCO population is dominated by the NS-NS binaries, 
causing at first a sharp rise followed by a slow decline at 
longer merger times (compare to the middle panel of Fig- 
ure 0). The fiattening of the distribution at 10^ - lO'* Myr 
is due to the small contribution of BH-NS and BH-BH sys- 
tems to the DCO population in this model. In the model 
shown in the middle panel (model F2), all dynamically 
stable MT episodes are assumed to be conservative. The 
distribution is rather fiat over a wide range of merger times 
and this is mainly a result of the assumption of conserva- 
tive MT and the wider orbits that form as a consequence. 
In the bottom panel, we show results from a model with 
low CE efficiency (model El). It is evident that the dis- 
tribution is depleted of short merger times. In contrast to 
the previous two models, the population is dominated by 
BH-NS and BH-BH systems, which have in general wider 
initial separations, and thus longer merger times than NS- 
NS binaries. Model H2 (no helium-star CE phases) shows 
a very similar behavior with the peak at short merger times 
disappearing. 

Other models lead to the double-peaked distributions 
similar to our standard-model results. The range of val- 
ues remains unchanged, but the relative strength of the 
two peaks vary slightly from model to model following the 
variations of the dominant classes within the DCO popu- 
lation. 



3.4.5. Double Compact Objects Merger Times 

Standard Model. The distributions of merger times (see 
eq. ^4|) of coalescing DCO are shown in Figure for our 
standard model. The double-peaked form is due to the two 
main populations of NS-NS systems and BH binaries and 
their double-peaked distributions of orbital separations. 
As already mentioned the NS-NS population is dominated 
by tight binaries formed through channels involving MT 
episodes from helium stars. These tight orbits imply very 
short merger times with a peak at ^ 0.3 Myr. Despite 
their higher masses BH-NS and BH-BH binaries are found 
to have merger times typical of ^ 1 Gyr, driven by their 
wider orbits and the stronger dependence of merger time 
on separation. Merger times also decrease with increas- 
ing initial eccentricities, leading to a stronger separation 
of the two peaks, since NS-NS binaries tend to have not 
only tighter orbits, but also higher eccentricities. 

In the bottom panel of Figure ^, we present the merger 
times of three different populations of coalescing NS-NS 
binaries (see §3.3.1): group I of tight, non-recycled NS- 
NS binaries; group II of tight binaries formed through the 
newly identified channels; and group HI of NS-NS formed 
through "classical" channels. It is evident that the large 
relative contribution of the tight NS-NS binaries drives 
the typical NS-NS merger times down to values close or 
below IMyr, whereas the "classical" NS-NS have merger 
times typical of the BH binaries. The identification of 
these short-lived binaries has important implications for 
the detectability of coalescing NS-NS binaries (see § 3.5.2). 

Both, Portegies Zwart & Yungelson (1998) and Fryer 
et al. (1999), presented merger time distributions of NS- 
NS binaries. They obtained typical times much longer 
than ours, close to 100-1000 Myr. This discrepancy origi- 
nates from the newly identified short-lived NS-NS systems. 
However, our classical population of NS-NS binaries have, 
as expected, merger times comparable to these found in 

^Following van der Kruit (1987) we have adopted a value of 2 X IO^^Lq for the blue luminosity of the Galaxy, the value used in several 
other studies, e.g., Phinney (1991) or Portegies Zwart & Yungelson (1998). Some recent results point to a value lower by a factor of ~ 2 (see 
discussion in Kalogera et al. 2001, and references therein). Had we adopted this lower value, all our coalescence rates would also be decreased 
by the same factor of ~ 2. 



3.5. 



3.5.1. 



Results on Double Compact Object Coalescence 
Rates 

Calactic Coalescence Rates of Double Compact 
Objects 



The calculated DCO coalescence rates have been cali- 
brated using the latest Type II SN empirical rates nor- 
malized to our Galaxy |^ (Cappellaro, Evans, & Turatto 
1999). In Table ^ coalescence rates are given for all the 
models in our study for each of the DCO classes as well 
as for the whole population. For our standard model (A), 
rates were calculated from a very large number of primor- 
dial binaries (3 x 10^) and sub-samples of DCO formed 
from each of 10^-binaries runs were used to examine and 
ensure the statistical accuracy of the results (see § 3.2). All 
other rates are based on models with iV > 10^ primordial 
binaries. For each model, we have also evolved an equal 
number of single stars, except for models Kl-2, for which 
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we appropriately adjusted the contribution of single stars 
based on the assumed binary fraction. 

In Figure |o[ we illustrate the rate dependence on the 
assumed kick velocity distribution. There is an overall 
decrease of rates with increasing kick velocity for every 
population of coalescing double compact objects, as the 
disruption probability for pre-coUapse systems increases 
with higher kick magnitudes. It is worth noting though 
that at small kicks (from zero to average magnitudes of 

30kms~^) rates are found to increase. This is due to the 
importance of kicks in creating tight binaries with merger 
times shorter than a Hubble time (see Fryer & Kalogera 
1997 and Fryer et al. 1999 in the NS-NS case). Rates re- 
main roughly constant for a — 20 ~ 50kms~^, as the two 
effects balance one another. It is also notable that the 
rate decrease is steeper for NS-NS binaries, and progres- 
sively flattens for BH-NS and BH-BH binaries. This is 
a result of the lower kick magnitudes assumed to be im- 
parted to BH. For very high kick magnitudes (Maxwellians 
with a > 400kms^^), BH-BH systems start dominating 
the population, as most of the NS-NS progenitors get dis- 
rupted. These results are in qualitative agreement with 
previous studies (e.g., Lipunov et al. 1997; Portegies Zwart 
& Yungelson 1998; Fryer et al. 1999). In Figure |l^, we 
also plot the rates from our standard model (kicks similar 
to the results by Cordes & Chernoff 1998) and the model 
B13 with "Paczynski-like" kicks. It is evident that, despite 
the very different shape of these distributions, they closely 
correspond to Maxwellian kicks with a ~ 240km s~^ and 
a ~ 150km s^^, respectively. This result once again in- 
dicates that the kick distributions are narrowly "filtered" 
by DCO binary properties: only binaries receiving kicks 
similar to their orbital velocities have a good SN survival 
chance (Kalogera 1996). The match between these two 
models and the specific cr-values of the Maxwellians in- 
dicates that the normalizations in the velocity range of 
interest happen to be very similar and the shape of the 
distributions outside this range becomes irrelevant. 

The effects of hyper-critical accretion combined with the 
assumed maximum NS mass can be understood on the 
basis of the results of models C and Dl-2. We find that 
hyper-critical accretion in CE phases not only increases 
the mass of the inspiraling compact objects and can con- 
vert NS to BH, but also leads to somewhat wider post-CE 
systems because part of the envelope is accreted and does 
not need to be expelled at the expense of orbital energy. In 
model C, we do not allow for any hyper-critical accretion 
and we find that the rate of NS-NS and BH-NS systems 
decreases whereas the rate BH-BH binaries remains un- 
affected. This combination indicates that it is the effect 
on the orbital period (increased rate of post-CE mergers) 
that dominates over the reduction of NS conversions to 
BH. BH-BH systems originate from wider progenitors and 
therefore are not much affected. The reduction in NS con- 
versions to BH turns out to be unimportant because the 
assumed maximum NS mass is rather high and the rate 
of conversions is low in the standard model (only 5% of 
all NS entering CE phases). In contrast, when the max- 
imum NS mass is reduced to 2M0 and 1.5 Mq (models 
Dl and D2), 34% and 80% respectively of NS entering CE 
phases collapse into BH. In these two models there is also 
a clear decrease in the rate of NS binaries and increase 



of the BH-BH coalescence rates, whereas the total DCO 
coalescence rates remain essentially constant (within our 
statistical accuracy, see §3.2). 

The effects of varying the effective CE efficiency (models 
El-3) are qualitatively similar to not allowing for hyper- 
critical accretion (model C). Coalescence rates tend to 
decrease with decreasing CE efficiency because of an in- 
creased rate of CE mergers. The main difference with 
model C is that the BH-BH rates are altered here. This is 
because varying the CE efficiency by factors of 2 or more 
affects post-CE binary separations much more than de- 
creasing envelope masses of an ejected envelope by a few 
tenths of a solar mass in model C. 

In models Fl-2 and Ll-2, we vary the parameters gov- 
erning mass loss in non-conservative MT episodes, and in 
particular in model F2, we obtain results for conserva- 
tive mass transfer. Most of the NS-NS progenitor sys- 
tems (~ 75%) start their mass transfer history with a 
non-conservative exchange, and therefore their coalescence 
rates might be altered. Also for many other compact bi- 
nary progenitors, non-conservative mass transfer phases 
may take place, particularly at the early stages of binary 
evolution, when the primary expands and evolves toward 
the red giant branch while its companion is still on the 
main sequence. However, non-conservative mass transfer 
phases do not drastically change the orbital separation, so 
we find that the overall coalescence rates do not change by 
much. Comparison to our standard model shows a deple- 
tion in coalescing BH-BH systems for both Fl-2 and Ll-2 
models, and also moderate decrease of NS-NS binaries for 
the Fl and L2 models. Coalescence rates of BH-NS bina- 
ries are not significantly altered. The similarities of the 
results come from the same effect that the fact that the F 
and L model parameters affect post-MT orbital separation 
in a similar way. In models F, we vary /„, the amount of 
material which is accreted onto the companion star, while 
in models L, we change the specific angular momentum 
j of the material lost from system during MT (1.0 — fa)- 
Both an increase of fa (F2) and decrease of j (LI) lead to 
an increase of the final post-MT orbital separations. On 
the other hand, both a decrease of fa (Fl) and increase of 
j (L2) lead to a decrease of the final post-MT separations, 
(see eq. || and eq. ^, or original study of Podsiadlowski et 
al. 1992)! Note, that orbital separations also depend on 
the mass ratios, but here we only discuss the changes rela- 
tive to the standard model. It turns out that these orbital 
separation changes have different effects on different DCO 
populations. For NS-NS progenitors, pre-MT orbits are 
already relatively tight, so in models Fl and L2 (F2 and 
LI) the NS-NS rate decreases (increases) as the frequency 
of mergers increases (decreases). Progenitors of BH-BH 
systems are more massive and generally have wider orbits 
at the onset of MT episodes compared to the progenitors 
of NS-NS binaries. Given the decreased post-MT separa- 
tions of model Fl and L2, the BH-BH rate drops as in the 
case of NS-NS systems and for the same reasons. How- 
ever, for the increased post-MT separations of model F2 
and LI, the BH-BH rate drops as well, this time due to the 
fact that many final BH-BH binaries are not tight enough 
to coalesce within the Hubble time. 

Variation of wind mass loss rates (models Gl and G2) 
affects more strongly the rates of BH binaries, which de- 
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crease with stronger mass loss. The extreme case of wind 
losses even stronger than in our standard model entirely 
eliminates the BH-BH population, but increases the over- 
all DCO rate, mainly because survival through CE phases 
and SN events is facilitated (lower-mass envelopes). 

As we pointed out in §2.1.2 (see also Belczynski & 
Kalogera 2001), the maximum mass of helium stars 
that develop convective envelope is somewhat uncertain 
(thought to lie within the range of 3.5 — 4.5 M©), and 
therefore we treat it as a model parameter. We find that 
reducing this value to Mconv = 4. DM© (model HI) re- 
sults in NS-NS coalescence rate decreased by 30%. This 
is expected since many of the NS-NS formation channels 
(see Table ^) involve CE evolution of helium stars, which 
is aborted in the absence of convective envelopes. If we 
further eliminate the possibility of CE evolution even for 
low-mass helium stars (model H2), the NS-NS coalescence 
rate drops to 0.9Myr~^. 

Although it is not thought to be relevant for our Galaxy, 
we examine one model where star formation is assumed to 
be burst-like and to have occurred lOGyr ago (model I). 
We find that all coalescence rates are increased because 
a larger fraction of them had enough time to evolve and 
have total lifetimes shorter than 10 Gyr. 

As expected, a flatter IMF (model J) favors the forma- 
tion of NS and even more of BH, so the contribution of 
BH-BH and BH-NS systems increases as the overall DCO 
coalescence rate increases as well. 

Varying the binary fraction affects the overall normaliza- 
tion of our population synthesis models. Binary systems 
mostly contribute to type Ib/c SN and therefore an in- 
creased binary fraction decreases the absolute number of 
Type II SN in the model and, for an assumed Type II SN 
normalization (based on the empirical estimates), leads to 
an increase of the coalescence rates. 

In model N, we see that rates of systems containing BH 
are similar to those of the standard model. However, the 
rates of NS-NS are much depleted. This is due to the fact, 
that in model N, with no helium star radial expansion, the 
NS-NS progenitors avoid the last CE episode of standard 
model formation channels. Therefore, many NS-NS sys- 
tems of model N are not tight enough to merge within the 
Hubble time. 

In models M, where we used different distributions for 
the initial binary mass ratios, the total DCO coalescence 
rates are reduced or increased, depending on whether ex- 
treme mass ratios (small q values) are favored or not. In 
model Ml the MT events are generally dynamically un- 
stable, leading very often to mergers long before compact 
object systems form. 

Evolution of model O influences only systems containing 
a BH, as in this model we have allowed for BH formation 
only through partial fall back. In the standard model, 
some BH were formed through a direct collapse of a mas- 
sive star, and no kick was imparted to such BH nor any 
material was lost from the binaries. Extending the influ- 
ence of BH formation through partial fall back increases 
the overall probability that binaries are disrupted. Thus 
coalescence rates of BH-NS and BH-BH binaries in model 
O are lower. 



3.5.2. Rate Increase Factors for Empirical Estimates of 
NS-NS Coalescence 

Coalescence rate estimates based on the observed sam- 
ple of close NS-NS systems need to account for all pos- 
sible observational biases acting against their detection. 
Systems with two NS are discovered as binary pulsars. 
The two coalescing systems that have been observed in 
the Galactic Observed sample of NS-NS systems contain 
recycled NS with long pulsar lifetimes (~ 10^ yr) and long 
merger times (> 10* yr), and empirical estimates are ob- 
tained under the assumption that the observed systems 
are representative of the Galactic population (Kalogera et 
al. 2001). 

As emphasized in our earlier work (Belczynski & 
Kalogera 2001), empirical coalescence rate estimates could 
be increased in light of the new short-lived NS-NS sub- 
populations identified in studies of their formation. Rate- 
increase factors can be calculated to account for the pop- 
ulations of (i) non-recycled NS-NS, and (u) tight NS-NS 
binaries with very short merger times formed through the 
new formation channels identified in the present study. Be- 
cause of their highly reduced lifetimes (pulsar and merger), 
the detection efficiency for these NS-NS sub-groups is prac- 
tically diminished (drops by ~ 3 orders of magnitude). 

In Table |, we present Galactic coalescence rates of: 
non-recycled NS-NS, tight NS-NS binaries with merger 
times shorter than 1 Myr, and the complete NS-NS popu- 
lation. In the last column of Table ^ we give the correction 
factors for empirical coalescence rates calculated from: 



Ti- — (^nr + 7?.tight) 

where TZ is the rate of all coalescing NS-NS binaries, TZm is 
the rate of coalescing non-recycled NS-NS (group I) , 7?.tight 
is the rate of systems (of group II and HI) with merger 
times smaller than 1 Myr. These correction factors imply 
an increase of the NS-NS rate estimated empirically based 
on the observed NS-NS with long pulsar and merger life- 
times. 

It is evident from Table ^ that systems formed through 
the new "helium-star" channels (group I and II) contribute 
the most to the correction factors. Our standard model 
prediction is that empirical rates may be increased by a 
factor of 2.5 to account for short-lived NS-NS pulsar sys- 
tems. Although the correction factors change significantly 
for a few of the models, i.e., from 1.4 to 6.1 for models 
with different CE efficiency (El-3), they remain roughly 
constant at ^ 2.5 for most of the models in our param- 
eter study. Only for the nonphysical model N, with no 
helium star radial evolution, the correction factor is found 
as expected to be 1.0. Note that even for model H2 (no 
helium-star CE evolution) there is a small fraction of short- 
lived systems (primarily due to kicks instead of CE orbital 
shrinkage) and the correction factor is 1.1. In model N no 
correction is needed, as in this model we do not form any 
NS-NS binaries in groups I nor II, and both these groups 
are primarily responsible for the increase factors (see also 
footnote in § 2.1.2). 

For our correction-factor estimates we have conserva- 
tively assumed (i) that a recycled NS is formed in all NS- 
NS binaries in which the first-born NS had an opportu- 
nity to interact with its non-degenerate companion, even 
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though this may not be true (given our hmited understand- 
ing of recychng we cannot be certain). We characterize 
non-recycled NS-NS only systems that evolved through 
a double CE phase involving two helium stars (channels 
NSNS:09, NSNS:11, and NSNS:12); (ii) that systems with 
merger times longer than 1 Myr can detected with the 
same efficiency as much longer lived binaries, even though 
in reality there is a continuum, given the continuous dis- 
tribution of merger times (Figure |^). Because of these 
two assumptions we should consider the derived upwards 
correction factors to be rather conservative and represent 
more of lower limits than rough estimates. 

3.5.3. Predicted Supernova and Star Formation Rates 

Cappellaro et al. (1999) estimated the rates of Type II 
SN and Type Ib/c SN to be 0.86±0.35 SNu and 0.14±0.07 
SNu for Sbc-d galaxies, where 1 SNu corresponds to one 
SN per 100 yr and the estimates are normalized to a blue 
luminosity of IO^^Lq^ . For an estimated Galactic blue 
luminosity of about Lq^ = 2 x IO^^Lq (van der Kruit 
1987; although sec footnote in § 3.5.1), we obtain 1.72 and 
0.28 SNu for Type II and Ib/c SN, respectively. The ratio 
of Type II to Type Ib/c SN turns out to be quite uncer- 
tain: 6.1 ± 4.0. Note that this ratio is independent of the 
assumed blue luminosity for our Galaxy. 

In our population synthesis calculations we keep track 
of all SN events. We use Type II events to normalize our 
models and then combine them with the Type Ib/c rates 
to derive the model-predicted Il-to-Ib/c ratios. We find 
them to lie in the range 1.8 — 3.7, for all the models in 
our parameter study, with the standard model yielding a 
ratio of 2.6. Although our predicted ratio values are con- 
sistent with the empirical estimates within the associated 
errors, we note that the theoretical values tend to be sys- 
tematically lower than the empirical estimates, implying 
that our Ib/c rates are rather high. In retrospect, this is 
actually expected given some of the assumptions in the 
MT events. We note that in dynamically stable and un- 
stable MT episodes, we assume that the donors are always 
stripped of their envelopes. However, it is not at all clear 
that this is a realistic assumption in small cases and it is 
quite possible that a small fraction of the H-rich envelope 
remains with the post-MT/CE donors. Since our classifi- 
cation of SN events as type Ib/c is based on the exploding 
stars having lost their H-rich and/or He- rich envelopes, it 
is possible that we are overestimating the Ib/c events be- 
cause of this assumption. We note, however, that all our 
models remain consistent with the empirical rates within 
the estimated errors. 

Another explanation for the possible overabundance of 
Type Ib/c SN in population synthesis model has been pro- 
vided by De Bonder & Vanbeveren (1998). Our results 
are consistent with their findings of Il-to-Ib/c ratios of 
2-3. The explanation they offered was related to possi- 
ble variations of massive binary formation in galaxies over 
a range of morphological types and to the empirical es- 
timates representing an average over many galaxies and 
not being appropriate for our Galaxy. However, in their 
study De Donder & Vanbeveren (1998) considered the re- 
sults obtained by (Cappellaro et al. 1993a, b) that were 
available at the time. In contrast, the most recent study 
of empirical SN rates (Cappellaro et al. 1999) accounts 



for different morphological types of galaxies, possibly im- 
plying that the MT assumptions (common to our and the 
De Donder & Vanbeveren 1998 study) may be primarily 
responsible for the calculated SN small rate ratios. 

Galactic SFR have been estimated to lie in somewhat 
broad ranges of 1 - 3MQyr-i (Bhtz 1997; Lacey & Fall 
1985) and - 1 - lOMgyr"! (Gilmore 2001). We use the 
calibration to the Type II SN rate by Cappellaro et al. 
(1999) with the adopted blue luminosity of our Galaxy: 
2 X 1O^°L0 (van der Kruit 1987) to calculate the Galac- 
tic SFR corresponding to our models. To do so we have 
to make an assumption about the extension of the IMF 
down to 0.08 Mq. If we assume that the IMF continues 
as a steep power-law down to the hydrogen-burning mass 
limit, we obtain values in the range 7- 28M0yr-i, which 
is clearly significantly higher than current estimates for 
the Galactic SFR. Based on the results of Scalo (1986, see 
also Kroupa, Tout, & Gilmore 1993), the IMF is thought 
to flatten for masses below about 1.0 Mq. With the as- 
sumed flattened IMF of Kroupa et al. (1993), we obtain 
much lower SFR values in the range of 3-9M0yr-i, with 
the standard model prediction of ~ 6M0yr~^. Moreover, 
had we adopted the lower, by a factor of two, blue luminos- 
ity of our Galaxy (see Kalogera et al. 2001, and references 
therein) the predicted SFR would decrease (by the factor 
of two) even further. We note the good agreement of the 
Galactic SFR estimates with our predictions for a flattened 
IMF. 

3.5.4. Comparison With Other Studies 

In § 3.4 we compared our results for the physical prop- 
erties of DCO populations to those of earlier studies wher- 
ever possible. In this subsection we explicitly focus on 
such a comparison based on results for coalescence rates 
from the following studies: Lipunov et al. (1997), Porte- 
gies Zwart & Yungelson (1998, hereafter PZY), De Donder 
& Vanbeveren (1998, hereafter DDV), Fryer et al. (1999, 
hereafter FWH), and Nelemans et al. 2001b. 

Lipunov et al. (1997) focused their study on the effect 
of SN kicks on coalescence rates of double compact ob- 
jects. Their calculations showed that, in general, the rates 
decrease approximately exponentially with increasing kick 
magnitude. This finding is in very good agreement with 
our calculations (see §3.5.1 and Figure We also see 
that the behavior of relative coalescence rates of NS-NS to 
BH-NS binaries is quite similar to what we obtain here. 
In both studies, the ratio of rates (NS-NS to BH-NS) de- 
creases with increasing average kick velocity. However, the 
absolute coalescence rates differ and this is easily under- 
stood as a result of different population synthesis assump- 
tions, important for the calculation of coalescence rates. In 
particular, Lipunov et al. (1997) use quite different com- 
pact object formation scenarios, as well as a different IMF. 

PZY studied the formation and evolution of NS-NS and 
BH-NS binaries and included a limited parameter study. 
Nelemans et al. (2001b) calculated populations of Galac- 
tic binaries with white dwarfs, NS, and BH, and estimated 
the low-frequency gravitational-wave emission of these bi- 
naries. They used a population synthesis code similar to 
that of PZY, with few modifications concerning BH for- 
mation. Although, many of the population synthesis as- 
sumptions differ, we find some similarities in the results. 



18 



In general, rates of NS-NS binaries are higher than those 
of the BH-NS systems. However, in all PZY models NS- 
NS dominate over BH-NS binaries, which is not true for 
all of our models mainly because of our more extensive pa- 
rameter study. Moreover, on average the rates obtained by 
PZY are smaller than ours, but for some comparable mod- 
els they differ only by factors of a few. This is due to the 
newly identified population of NS-NS binaries, as well as 
to the other differences in model assumptions. PZY note 
that due to hyper-critical accretion onto NS in CE phases, 
the population of BH-NS may dominate over NS-NS sys- 
tems. Although, they do not perform actual calculations 
of hyper-critical accretion events (they assume that NS al- 
ways collapse to BH), their results are in agreement with 
ours, provided that the maximum NS mass is rather low 
(smaller than 2.0 M©, see rates of models Dl and D2 in 
Table ^) . The predicted BH-BH coalescence rates are di- 
minishingly small in both studies (BH-BH binaries form 
but in very wide orbits with merger times longer than a 
Hubble time). This is the combined result of a number 
of factors: in these two studies BH are assumed to form 
(i) symmetrically and no birth kicks are taken into ac- 
count, and (ii) only from progenitor stars more massive 
than 40 Mq , which biases the population to lower frequen- 
cies and wider orbits. 

With respect to the DDV study, we note the difference 
between formation rates that correspond to the entire pop- 
ulation and coalescence rates that correspond to the sub- 
group with coalescence times shorter than a Hubble time. 
DDV reported formation rates and hence a direct com- 
parison with their results is difficult. Overall their forma- 
tion rates are very high (not only compared to our results 
but also all the other population studies): ~ 400 NS-NS, 
~ 2000 BH-NS and ~ 140000 BH-BH systems per Myr in 
Galaxy. DDV comment that their BH-BH formation rate 
is surprisingly high. One reason may be their calibration 
method based on an assumed massive star formation rate 
of one massive star per year. This does not easily trans- 
form into a star formation rate as used in this work and 
consequently makes the comparison very difficult. 

FHW focused their work to the study of gamma-ray 
burst progenitors and among them NS-NS and BH-NS 
binaries. They performed quite an extensive parameter 
study and calculated coalescence rates for all their mod- 
els. Once again, we note the very similar dependence of 
their rates and ours on the average kick magnitude. One 
striking difference comes with the very wide ranges of their 
predicted rates (e.g., their 4 orders of magnitude change 
of NS-NS coalescence rate compared to about 2 orders of 
magnitude found here). The enormous spread of rates in 
the FHW work comes from the fact, that in one of their 
models they, non-physically, increased (by factor of four) 
the maximum stellar radii to explore the uncertainty re- 
lated to radius determination in evolved stars. Such a 
change significantly affects the evolution of many binaries, 
since the stellar radius is a crucial quantity in judging on 
the occurrence or non-occurrence of one of the most impor- 
tant binary processes: MT events. The radii we calculate 
for stars come from the stellar evolution models, and they 
depend on the star mass, its evolutionary stage, and its 
composition, and they are not altered in any nonphysical 
way. However, the FWH models of varied radii show the 



sensitivity of the results of population synthesis to this 
parameter, which is not very well constrained by stellar 
evolution models. 

The results on coalescing compact object binary rates 
are in good qualitative agreement with previous theoreti- 
cal predictions. Although we have noted some systematic 
differences we are able to attribute them to the different 
population synthesis model assumptions, as well as to the 
effects of the newly identified populations of NS-NS bina- 
ries. 

4. DISCUSSION 

Close binaries of NS and BH have attracted an increased 
interest in recent years primarily because of their con- 
nection to gravitational-wave detection and possibly to 
gamma-ray burst progenitors. Here we focus the discus- 
sion of our results in the context of gravitational-wave 
detection by the upcoming ground-based interferometers 
(e.g., LIGO) and the prospects of detecting inspiraling 
DCO. An analysis of these and other populations in con- 
nection to gamma-ray bursts is presented in Bclczynski, 
Bulik & Rudak 2002c and Pcrna & Belczynski 2002. 

Our motivation in initiating this study was to examine 
DCO populations in view of some developments in the un- 
derstanding of CE evolution, particularly the possibility of 
double CE ejection (suggested by Brown 1995) and hyper- 
critical accretion (suggested and quantified by Chevalier 
1989, 1993; Brown 1995; Bethe & Brown 1998, and in the 
present study), and the evolution of low- mass helium stars. 
Some, although not all, of these possibilities have been ex- 
plored in some of the earlier population synthesis studies, 
but not in a consistent and equally detailed way (mostly 
some extreme cases have been examined; for example, ei- 
ther hyper-critical accretion has been ignored or all NS 
entering a CE have been assumed to collapse to BH) . Our 
goal is primarily two-fold: (i) to examine the predicted 
rates for various DCO classes, focusing mainly on the rel- 
ative formation frequencies and their behavior with a large 
number of model parameters, and (ii) to examine in de- 
tail the physical properties of the various populations, the 
origin of their characteristics, and the links to certain key 
evolutionary effects or phases, and to identify the most 
robust of the qualitative features in the distributions of 
binary properties. 

In the course of our investigation we came across a num- 
ber of new formation channels leading to rather efficient 
formation of NS-NS binaries. Apart from an increase in 
the predicted rates, examination of their properties re- 
vealed that in their majority these NS-NS binaries form 
a very distinctive class with tight orbits and short merger 
timescales (i.e., lifetimes). The common thread connect- 
ing all these new formation paths is the evolution of low- 
mass helium stars and the implications of the fact that 
they can develop partially or fully convective envelopes. 
We find that a large number of MT episodes initiated by 
these evolved stars lead to high survival rate through a 
CE phase and the production of very tight binaries, which 
in turn have very high survival probability through the 
second core-collapse event. Further, following the basic 
argument made by Brown (1995) we have allowed for the 
possibility of two giant stars (He-rich more relevant to ef- 
ficient DCO formation, but also H-rich) involved in the 
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dynamically unstable MT to eject the two combined com- 
mon envelopes. It is these possibilities that lead to qualita- 
tive changes in the DCO population characteristics, which 
have been discussed above. 

It is important to acknowledge that the viability of these 
new formation channels has not been examined with de- 
tailed evolutionary and hydrodynamical calculations. De- 
spite the difficulties of fully understanding the details of 
CE phases, questions such as: "are DCE events realistic?" , 
"can evolved helium stars survive a CE phase?" , and "un- 
der what conditions?" , are currently traceable (model 112 
corresponds to the case of the answers to the above ques- 
tions being no for helium stars). It seems that such in- 
vestigations are necessary before we can include these new 
channels as part of what is considered "standard" ways of 
forming DCO. Nevertheless it is interesting to examine the 
implications of such evolutionary phases (see also Fryer et 
al. 1999; Nelemans et al. 2001a). 

In the context of gravitational-wave detection from in- 
spiraling NS-NS binaries, these relatively short-lived sys- 
tems imply that empirical rates derived based on the ob- 
served sample (with much longer lifetimes) should possibly 
be raised by factors of 2-3 typically. The reason is that 
they correspond to a Galactic NS-NS population that is 
not represented in the observed binary pulsar sample but 
would very well contribute to inspiral events. Using the 
results of Kalogera et al. (2001) for the empirical NS-NS 
coalescence rate, we find that their most optimistic pre- 
diction for the LIGO I detection rate could be raised to at 
least 1 event per 2 years, and their most pessimistic LIGO 
II detection rate could be raised to 5 events per year or 
even higher. 

In Table ^ and Table |^ we present our theoretical pre- 
dictions of detection rates of different binary merger types 
for LIGO I and LIGO II, respectively. These detection 
rates correspond to our coalescence rates of double com- 
pact objects calculated with the StarTrack population syn- 
thesis program. Using the extragalactic extrapolation of 
Kalogera et al. (2001) and the maximum sensitivity dis- 
tances of LIGO I and LIGO II for a given binary merger 
type (see Kalogera & Belczynski 2001), we converted our 
Galactic coalescence rates to detection rates. We find that 
the LIGO I detection rates are quite low with a maximum 
total rate of a couple events per year and with rates signif- 
icantly below 1 event per year for many models. However, 
we find that the prospects for double compact object inspi- 
ral detection are very encouraging for LIGO II, for which 
we predict at least ~ 10 detections per year, even in the 
most pessimistic case. Moreover, for LIGO II, the total 
detection rate of compact object binaries is as high as few 
hundred events per year for most of our models. 

The uncertainties of population synthesis method (re- 
flected on the predicted detection ranges) seem much re- 
duced compared to the earlier results (see Fryer et al. 1999; 
Kalogera et al. 2001). This is related to the fact that in 
the present study we adopted a given set of single star 
evolution models (however note that we did explore the 
effects of varying wind mass-loss rates and helium-star 
evolution), and we have adopted a physically motivated 
initial-remnant mass relation based on hydrodynamical 
calculations of core-collapse events. In light of some in- 
dicative results from rotating star models (Heger, Langer 



& Woosley 2000; Heger & Langer 2000), we can expect 
stellar evolution models to be updated in the future and 
we would regard this as progress in the field. 

We regard the investigation of the physical properties of 
DCO as important as that of the rates and quite revealing 
in terms of their origin and robustness. The details have 
already been discussed earlier in the paper. From the point 
of view of gravitational- wave astrophysics, we hope these 
results will open new directions in anticipating and under- 
standing properties of gravitational-wave sources. For ex- 
ample, distributions over components masses can be used 
to produce fake inspiral data to test the data analysis tools 
currently under development or evaluate detection efficien- 
cies (we are currently involved in such an activity already 
planned within the LIGO Scientific Collaboration as part 
of the preparation for the LIGO I Scientific run in 2002). 
Predictions for physical properties can be used in develop- 
ing specialized data-analysis tools to explore an astrophys- 
ically motivated parameter space, in cases where the unre- 
stricted parameter space of inspiral signals is just far be- 
yond any current or near-future computational capabilities 
(see Kalogera 2000). Looking a little more into the future, 
in the era when gravitational-wave astronomy is possible 
and physical properties are measurable, the identification 
of strong qualitative features (such as those discussed here: 
relative contributions of DCO classes, or presence or ab- 
sence of certain peaks and other features) will allow us to 
evaluate some of the surviving theoretical models. 
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APPENDIX 

HYPER-CRITICAL ACCRETION ONTO COMPACT OBJECTS DURING COMMON ENVELOPE PHASE 

Let us denote the mass of the eompaet object (NS or BH) by Ma, the mass of its companion (H-rich or He-rich giant) 
by M-Q, its core mass by Mb. core it^ radius by Rb.cc: f^-nd the binary semi- major axis by A^c- Accretion onto the compact 
object can be initiated only once the binary separation becomes equal to the radius of the expanding giant donor. Both 
quantities have changed since the time of Roche-lobe filling to ^acc and -Re.acc) so that ^acc = -RB,acc- CE evolution and 
accretion onto the compact object will end when the giant's envelope has been ejected and the donor's mass changes from 
Mb,! to MB,f = Mb, core- Throughout the derivation ace is the CE efficiency parameter and A is the numerical factor 
scaling the binding energy of the donor. The orbit is expected to be circular due to circularization on very short time 
scale as the donor approached Roche-lobe filling. 

Following Bethe & Brown (1998) we write the energy loss rate related to the accretion onto the compact object 

E^cc = 0.5cdG(MB + Ma)A-^Ma (A1) 

where Cd is the drag coefficient of the compact object with respect to the donor's envelope, and it is assumed Cd = 6 
(Shima et al. 1985). The rate of orbital energy dissipation is given by: 

- Eorh = 0.5GMbA-^Ma + 0.5GMaA-^Mb - 0.5GMaMbA-'^A (A2) 

Orbital energy is dissipated due to the dynamical friction of accreting neutron star in the giant envelope. Thus, we 
may compare equations A.l and A. 2 (.Bacc = — -^orb) and after taking derivative in respect to dMe we obtain, 

[cd(MB + Ma) - Mb] ^ = -MaMb + Ma (A3) 

During the phase of accretion we may express binding energy of donor envelope as, 

- Eund = GMb(Mb - MB,coro)A-M-i (A4) 

where instead of the donor radius we use the binary separation, which at the start of the accretion phase is equal to the 
donor radius (^acc = -RB.acc)- The rate of binding energy change can be written as, 

- £;bi„d = GX-^A-\2Mb - MB,corc)MB + GMBA-lA-2(MB,core " Mb) A (A5) 

The donor envelope is ejected on the expense of the binary orbital energy, with an efficiency described by parameter 
ace- Thus, we write CE energy balance with respect to the donor mass as, 

-"-(-dMB^ = ^"^^ ^""'^ 

We use equations A. 2, A. 5, and A. 6 to write: 

= MBA-i[-2A-ia-i(MB,core " Mb) + Ma]-^ - 2A-ia-i(2MB - Mg.core) - Ma (A7) 

We can next write out two ordinary differential equations to be solved: 

dA _ MaMbK^ + 2X-'^a-J'{2MB-MB,coTe) + MA 



dMB MaMIA-^H^^ + MbA-1 [-2A-iac-e'(MB,core - Mb) + Ma] 
dMA 



(A8) 



dMB 



MaK\1 - MbA- Vi(Ma, Mb, ^)) (A9) 



where h\ = Cd (Mb + Ma) — Mb , and /i (Ma , Mb , A) is equal to the right hand side of equation (A. 8) . We know the initial 
and final values for Mb, the initial value for Ma, and we can calculate the initial value for A when accretion is initiated. 
To do so we use the CE energy balance for the pre-accretion period (— aceAi^orb = AE'bind) and obtain: 

^ ^ 2(Mb - MB,core) + AaceMA 

2(Mb - MB,core)i?i' + AaceMA^cV ^ ' 

where we set i?B.acc = ^acc- Equations A. 8 and A. 9 describe the accretion phase during CE evolution [A < Aacc), and 
we integrate them numerically from Mb,! to Mb, core to obtain the final binary separation and final mass of the accreting 
compact object. 
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Table 1 

Population Synthesis Model Assumptions 



Model Description 



A 




standard model described in § 2.1 and § 2.2 


Bl- 


-13 


zero kicks, single Maxwellian with 

a = 10,20,30,40,50, 100, 200, 300, 400, 500, 600 kms-\ 

"Paczynski" kicks with a = 600kms~^ 


C 




no hyper-critical accretion onto NS/BH in CEs 


Dl- 


-2 


maximum NS mass: Mmax.NS = 2, 1.5 Mq 


El- 


-3 


acE X A = 0.1,0.5,2 


Fl- 


-2 


mass fraction accreted: fa = 0.1, 1 


Gl- 


-2 


wind changed by /wind = 0.5,2 


Hl- 


-2 


Convective Helium giants: Afconv = 4.0, Mq 


I 




burst-like star formation history 


J 




primary mass: oc Mj~^'^^ 


Kl- 


-2 


binary fraction: /bi = 0.25,075 


Ll- 


-2 


angular momentum of material lost in MT: j = 0.5, 2.0 


Ml 


-2 


initial mass ratio distribution: oc q~^'^ 


N 




no heUum giant radial evolution 







partial fall back for 5.0 < Mqo < 14.0 Mq 



Table 2 

Population Synthesis Accuracy 



Coalescing 


Mean 


Max Change 


Mean 


Population 


Rate'^ 


in Rate** 


Number 


NS-NS 


52.7 ± 1.1 (2%) 


4.4 (8%) 


1754 


BH-NS 


8.1 ± 0.6 (7%) 


2.2 (27%) 


269 


BH-BH 


25.6 ± 0.8 (3%) 


3.9 (15%) 


852 


Total 


86.3 ± 1.2 (1%) 


5.5 (6%) 


2875 



'^Galactic Coalescence Rate (Myr ^). 
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Table 3 

Double Compact Object Formation Channels - Standard Model 



Evolutionary History 



Formation 


Relative 




njinciciicy 


NSNS:01 


20.3 % 


NSNS:02 


10.8 % 


NSNS:03 


5.5 % 


NSNS:04 


4.0 % 


NSNS:05 


3.2 % 


NSNS:06 


2.5 % 


NSNS:07 


2.2 % 


NSNS:08 


2.0 % 


NSNS:09 


2.0 % 


NSNS:10 


1.6 % 


NSNS:11 


1.5 % 


NSNS:12 


1.5 % 


NSNS:13 


1.0 % 


NSNS:14 


3.0 % 


BHNS:01 


4.5 % 


BHNS:02 


1.6 % 


BHNS:03 


1.3 % 


BHNS:04 


2.0 % 


BHBHrOl 


17.7 % 


BHBH:02 


10.5 % 


BHBH:03 


1.4% 



NC:a^b, SN:a, HCE:b^a, HCE:b^a, SN:b 
NC:a^b, SCE:b^a, NC:a^b, SN:a, HCE:b^a, SN:b 

SCE:a^b, SN:a, HCE:b^a, HCE:b^a, SN:b 
NC:a^b, SCE:b^a, SCE:b^a, SN:b, HCE:a^b, SN:a 

DCE:a^b, SCE:a^b, SN:a, HCE:b^a, SN:b 
SCE:a^b, SCE:b^a, NC:a^b, SN:a, HCE:b^a, SN:b 
NC:a^b, NC:a^b, SN:a, HCE:b^a, HCE:b^a, SN:b 
NC:a^b, DCErb^a, SN:a, HCE:b^a, SN:b 

DCE:a^b, DCE:a^b, SN:a, SN:b 
NC:a-*b, SCE:b-^a, SN:b, HCE:a^b, SN:a 
NC:a^b, SCE:b^a, DCE:b^a, SN:a, SN:b 
NC:a^b, SCErb^a, DCEra^b, SN:a, SN:b 
DCE:a^b, SN:a, HCE:b^a, SN:b 
all other 

NC:a^b, SN:a, HCE:b^a, SN:b 
NC:a^b, SCE:b^a, SN:a, SN:b 
SCE:a^b, SN:a, HCE:b^a, NC:b^a, SN:b 
all other 

NC:a^b, SN:a, HCE:b^a, SN:b 
NC:a^b, SCE:b^a, SN:a, SN:b 
all other 



'^Normalized to the total DCO population. 

''Sequences of different evolutionary phases for the primary (a) and the secondary 
(b): non-conservative MT (NC), single common envelope (SCE), double common 
envelope (DCE), common envelope with hyper-critical accretion (HCE), supernova 
explosion/core-collapse event (SN). Arrows mark direction of MT episodes. 
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Table 4 

Galactic Double Compact Object Coalescence Rates 



ModeP 


NS-NS 


BH-NS 


BH-BH 


Total 


A 


52.7 


8.1 


25.6 


86.3 


Bl 


292.4 


18.2 


32.7 


343.2 


B2 


299.6 


19.4 


31.8 


350.8 


B3 


302.2 


19.6 


34.2 


356.0 


B4 


285.2 


19.1 


34.2 


338.5 


B5 


251.0 


19.5 


34.3 


304.7 


B6 


226.8 


16.4 


34.1 


277.3 


B7 


128.1 


14.6 


30.7 


173.3 


B8 


57.5 


10.1 


29.2 


96.8 


B9 


33.2 


5.7 


23.2 


62.1 


BIO 


18.2 


3.7 


21.0 


42.9 


Bll 


12.0 


2.1 


18.1 


32.2 


B12 


8.0 


1.6 


15.1 


24.6 


B13 


91.0 


10.3 


27.3 


128.6 


C 


43.2 


5.6 


23.2 


72.1 


Dl 


33.6 


23.3 


31.1 


88.0 


D2 


9.1 


36.2 


42.2 


87.5 


El 


2.7 


4.8 


5.6 


13.1 


E2 


23.5 


6.3 


23.1 


53.0 


E3 


109.0 


8.7 


11.5 


129.2 


Fl 


22.1 


9.3 


8.7 


40.1 


F2 


54.3 


8.6 


7.2 


70.1 


Gl 


43.9 


14.2 


75.6 


133.7 


G2 


92.2 


1.3 


0.0 


93.5 


HI 


37.9 


7.8 


26.6 


72.3 


H2 


0.9 


6.0 


26.3 


33.2 


I 


54.5 


10.0 


33.6 


98.1 


J 


58.1 


12.8 


41.9 


112.8 


Kl 


22.5 


3.4 


10.4 


36.2 


K2 


90.2 


13.5 


41.6 


145.4 


LI 


78.9 


9.2 


10.0 


98.1 


L2 


12.0 


6.2 


10.5 


28.7 


Ml 


6.2 


4.0 


5.8 


16.0 


M2 


114.2 


8.4 


31.5 


154.1 


N 


34.4 


10.7 


24.5 


69.6 





51.9 


5.7 


4.0 


61.6 



for definition of models see Table 



Table 5 

NS-NS Empirical Coalescence Rate Correction Factors 





„ T 


( Tj^vi 1 r~i 1 1 
vjriUUjJ ii 


Prmin TTT 


±ULcll 


JZL/llipii ICcLl XXclLtl 




till 


Tmcrg ^ ^ iviyr 


Tmorg ^ ^ iViyr 


XVcLLc 




A 


4.4 


26.9 


0.1 


52.7 


2.5 


Bl 


6.0 


146.7 


0.0 


292.4 


2.1 


B6 


6.5 


121.6 


0.0 


226.8 


2.3 


B7 


5.2 


68.0 


0.1 


128.1 


2.3 


B8 


4.7 


27.6 


0.2 


57.5 


2.3 


B9 


4.3 


17.0 


0.1 


33.2 


2.8 


BIO 


2.8 


8.8 


0.1 


18.2 


2.8 


Bll 


2.2 


6.1 


0.0 


12.0 


3.3 


B12 


1.9 


4.2 


0.0 


8.0 


4.2 


B13 


4.6 


48.4 


0.1 


91.0 


2.4 


C 


3.2 


22.1 


0.2 


43.2 


2.4 


Dl 


4.9 


15.1 


0.0 


33.6 


2.5 


D2 


3.6 


2.8 


0.0 


9.1 


3.3 


El 


0.4 


1.9 


0.0 


2.5 


6.1 


E2 


3.1 


15.1 


0.1 


23.5 


4.5 


E3 


5.2 


28.0 


0.2 


109.0 


1.4 


Fl 


2.3 


12.2 


0.1 


22.1 


2.9 


F2 


2.3 


16.1 


0.8 


54.3 


1.5 


Gl 


3.3 


23.7 


0.1 


43.9 


2.6 


G2 


7.3 


38.0 


0.2 


92.2 


2.0 


HI 


3.3 


20.9 


0.0 


37.9 


2.8 


H2 


0.0 


0.0 


0.1 


0.9 


1.1 


I 


4.0 


28.5 


0.0 


54.5 


2.5 


J 


4.4 


29.4 


0.1 


58.1 


2.4 


Kl 


1.8 


11.2 


0.1 


22.5 


2.5 


K2 


7.4 


44.9 


0.2 


90.2 


2.5 


LI 


6.3 


32.7 


1.1 


78.9 


2.0 


L2 


2.0 


5.6 


0.1 


12.0 


2.8 


Ml 


0.2 


3.9 


0.0 


6.2 


3.0 


M2 


14.0 


51.6 


0.2 


114.2 


2.4 


N 


0.0 


0.0 


1.1 


34.4 


1.0 





4.3 


26.5 


0.1 


51.9 


2.5 



''for definition of models see Table 



Table 6 

Predicted LIGO I Detection Rates (yr^^) 



Binary 


Standard 


Range 


Type 


Model 


(all models) 


NS-NS 


1 X 10-2 


2 X 10-4 - 7 X 10-1 


BH-NS 


2 X 10-^ 


2 X 10-3 - 7 X 10-2 


BH-BH 


8 X 10-^ 


0-2 


Total 


8 X 10-1 


2 X 10-3 - 2 
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Table 7 

Predicted LIGO II Detection Rates (yr~^) 



Binary 


Standard 


Range 


Type 


Model 


(all models) 


NS-NS 


6 X 10^ 


1 - 4 X 10^ 


BH-NS 


8 X IQi 


9 - 4 X 10^ 


BH-BH 


2 X 10^ 


- 8 X 10^ 


Total 


3 X 10^ 


1 X 10^ - 8 X 10^ 



No FB 



I 1 15 

O 



Partial ¥B Partial 1<'B - 

Complete FB Complete FB ; 




15 7 




'a 




15 7 

10 7 

5 7 
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Fig. 1.-- Final remnant masses shown as a function of progenitor ZAMS mass, for single Population I stars (Z = 0.02) 
and for three different wind mass-loss rates. Top panel: standard HPT wind mass-loss rate. Middle panel: Wolf-Rayet 
wind rates decreased by two. Bottom panel: all wind mass-loss rates decreased by a factor of two. Different symbols 
represent different remnant types (square: white dwarf; cross: neutron star; triangle: black hole). Stellar masses just prior 
to the collapse are shown by the dot-dashed line. In the top panel, we also mark the ranges of initial masses, for which 
NS or BH are formed without any fall back, and with partial or complete fall back in core-collapse evcnits of massive stars. 
In all panels, a maximum NS mass of 3 Mq has been assumed. Masses of WD remnants are calculated as in HPT. 
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Fig. 2. — Normalized distributions of the eompact object masses in coalescing NS-NS, BH-NS and BH-BH systems for 
the standard model. We show the mass distributions for first- and second-born compact objects in the top and bottom 
panel, respectively. The distributions are calculated using a total number of 28655 coalescing DCO formed out of 10^ 
primordial binaries. Note that the standard model choice of the maximum NS mass is: M^f^ = 3 Mq (marked with the 
dashed line). 
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Fig. 3. — Normalized distributions of the compact object masses in coalescing NS-NS, BH-NS and BH-BH systems for 
model El with very low CE efficiency: acE x A = 0.1. We show the mass distributions for first- and second-born compact 
objects in the top and bottom panel, respectively. The distributions are calculated using a total number of 1283 coalescing 
DCO formed out of 3 x 10^ primordial binaries. This case of very low CE efficiency results in very low DCO rates and 
relatively small-number statistics. Note that the standard model choice of the maximum NS mass is: M^-^ = 2>Mq 
(marked with the dashed line). 
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A [Ro] 

Fig. 4. Normalized distributions of orbital separations of eoalescing DCO at the time of their formation, for our standard 
model, for all DCO (top) and for the three classes separately (bottom). For comparison, we also show the observed orbital 
separations of the 2 coalescing NS-NS systems found in the Galactic field, PSR 1913+16 and PSR 1534+12. 




Fig. 5. — Normalized eccentricity distributions of coalescing DCO for our standard model at the time of their formation 
(dot-dashed line) and at the time the orbit decayed so that the corresponding gravitational- wave frequency is about 40 Hz, 
i.e., the system is entering the LIGO I band (solid line). For comparison wc also show the two observed short-period 
NS-NS systems: their current (open symbols) eccentricities and the eccentricities as they enter the LIGO I band (filled 
symbols). 
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Fig. 6. — Normalized distributions of DCO systemic velocities acquired after the first core-collapse event, for our standard 
model, for the whole population (top) and for each DCO class separately (bottom). 
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Fig. 7. — Normalized distributions of DCO systemic velocities acquired after the second core-collapse event, for our 
standard model, for the whole population (top) and for each DCO class separately (bottom). 
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Fig. 8. — Normalized merger time distributions shown for our standard model, for the whole DCO population (top), 
for each of the DCO classes (middle), and for each of the three major NS-NS groups (bottom). The merger of the two 
observed systems are also shown in the top panel. 
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Fig. 9. — Normalized merger time distributions of coalescing DCO for models Bl (top panel), F2 (middle panel) and El 
(bottom panel). The distributions are normalized to the total number of coalescing DCO: 10841, 2352, 1283 for models 
Bl, F2 and El, respectively; formed out of 10^ primordial binaries for model Bl and F2, and out of 3 x 10^ binaries for 
model El. 
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Fig. 10. — Dependence of Galactic coalescence rates on the assumed natal kick velocity distribution. Lines connect rates 
for models B1-B12 and the horizontal scale shows the width of Maxwellian kick distribution of a given model. Points 
mark rates for our standard model (A) and model with "Paczynski-like" kick distribution (B13). The horizontal position 
for these two models does not correspond to the horizontal axis scale, but instead is chosen so that the predicted rates 
approximately match the curve for the Maxwellian kick distributions. 



